//  scp69q.cpp Dec 8  looking at CdS in GeO2  Tanaka...
//  69p Dec 8 (L=0 only)
//   [ goes with 67p and 68p (toroidal) ]
//  69  Dec  L=0 version
//     Amplitude of B is also plotted as check of energy conservation.
//  68  Dec 6 toroidal version:
//   tested for L=1 and L=2 for light matrix
//  As an additional check note that B has a constant amplitude of 1,
//  which follows from energy conservation.
//   - I notice there is a "rotational" peak for Si in SiO2.  This
//     may explain the downshift of the depolarized spectrum.
//  67  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 object 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
//  59  looking for bug
//  58  relatively strong peaks before the beta=2.64 peak.
//  57  anomalous peak at beta=3.0 for these values (again) A & B
//  56  tried double rather than double.
//    Floating precision only seems to matter at beta below 1, where
//    there is slight jaggedness of the curve.
//    - plot both A and B.
//    - found bug in evalobj() - corrected error
//  55  This one shows structure even when matrix equals particle.
//       I don't know why this structure exists.
//   54  This one gives anomalous peak at beta=2.99 !! ( i don't know why)
//   53  L=2 spheroidal transverse incident wave
//   This gives exact values including 4.866 this time.
//   It is also interesting to note the several small peaks that
//   come before the peak at beta=2.64.  I am not sure what causes these.
//   52    exchanged A and B  (also matches exact values)
//  51  try L=2   This version does seem to give correct peaks
//      at exact values:  beta=2.64 [not 4.866] 8.330 9.781
//      It started to work correctly after I exchanged var's 6 and 7
//      in the linear solution procedure.
//  50   I notice a broad lump at low beta which I believe is due to
//      the linear translational mode of vibration.
//  49  L=2 spheroidal longitudinal (A only -- should also show B)
//   - Does closely match exact results of beta=3.425, 6.772, 7.746, 10.696
//  48
//  47  check if aux var's in use, and release when done.
//  46  L=1 toroidal, works for L=2 as well  (checked against exact beta)
//    L=1:  5.8, 9.1
//   - add factor of beta to account for phonon density of states.
//  45  equal velocity, equal density case -- no peak, as expected
//   Also test of correctness for nu=0.25 (exact = 4.44, 10.494)
// scp44 Dec 4  functional; makes plot of A versus beta for L=0
//    (corrected error where k factors were omitted)  R invar. check
// scp41  Dec 4  stress(,,,int imater)
//  40  Dec 4   try L=0 problem...
//  39  PARTIC, MATRIX, LONGIT, TRANSV
//  38
//  37  try to move cos() up by 1    maybe working...
//   simplify dispterm()
//  36  sphbessel()
//  35   NMATER = number of materials
// scp34 Dec 3  rearranged to allow generalization to multi layers
// scp33 Dec 1 improved dispterm(); seems to work
//  31
// scp30  Dec 1  ingoingwave outgoingwave()  multbyA()... subtract()
//  Just now it occurred to me that there are two different kLp's and kTp's!!
//  There is a kLp inside and outside the sphere  kLpp and kLpm!!!
//  So I need to re-engineer this again!
//  29   L=2 toroidal
//  28   L=2 transverse;   settozero();  multby() zero case
//  [note 27 got corupted by editor]
// scp26  L=2 longitudinal case  ur, uth and uph
//  25  stress()
//  24   fixed error in strain(); error check in add()
//  23  Nov 29  multby() multbyr() ...
//  22  calculate six components of strain
// scp21  test: divergence of curl of arbitrary vector field is zero.
// scp20  curl
//  19
// scp18 Nov 28 demo of zero laplacian; extend legen to order 5
//  17 works
//  16  minor clean up, generalize NFAC
//  scp15  legendre, sph. bessel of 1st and 2nd kind
//  scp14  Nov 28  multiply()
// scp13 works
// scp12 works
//  11 works
//   9
// scp8   cleanup function
// scp7   gradient, divergence (allowing Laplacian of scalar field)
//   (but there is no simplification of terms)
// scp6   Nov 27   theta derivative too. (minor bugs fixed)
// scp5  j(2)(kr) P(2)(costheta)  (order 2 or less only so far)
//   ... so far this can handle scalar fields.  It might be
//   useful to extend it to a general vector field.
// scp4  sph. Bessel function of order L   Mon Nov 25
// 3 // calculate Bessel function of order 2
//  2 // r derivative ddr(,)
// scp1.cpp  Sin Cos Power - for Sph. Bess. func.  Nov 25, 2002
// 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

// nm[0][]   sin(theta)
// nm[1][]   cos(theta)
// nm[2][]   sin(r)
// nm[3][]   r
// nm[4][]   cos(r)
// nm[5][]   kLp (longitudinal speed of sound in particle)
// nm[6][]   sin(kLp r)
// nm[7][]   A
// nm[8][]   cos(kLp r)
// nm[9][]   kTp (transverse speed of sound in particle)
// nm[10][]  sin(kTp r)
// nm[11][]  B
// nm[12][]  cos(kTp r)
// nm[13][]  kLm (longitudinal speed of sound in matrix)
// nm[14][]  sin(kLm r)
// nm[15][]  C
// nm[16][]  cos(kLm r)
// nm[17][]  kTm (transverse speed of sound in matrix)
// nm[18][]  sin(kTm r)
// nm[19][]  D
// nm[20][]  cos(kTm r)
// . . .
// nm[NFAC-1][]  ...

#include <iostream.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>
#include <dos.h>
#include <graphics.h>
#include <complex.h>
#define NOBJ 19
#define NTRM 3600
#define NMATER 2
#define NFAC 5+8*NMATER

#define ntrm0 1200
#define ntrm1 1200
#define ntrm2 1200

double far cr[NTRM],ci[NTRM];
int far nm0[NFAC][ntrm0];
int far nm1[NFAC][ntrm1];
int far nm2[NFAC][ntrm2];

int far bt[NOBJ+3]; // base term
int far nt[NOBJ+3]; // number of terms
int far mnt[NOBJ+3]; // max. number of terms
void dispterm(int it2,int ic,int ir);
void ddr(int iobj2,int iobj1); //  Take Deriv of iobj1 with respect to r
void ddth(int iobj2,int iobj1); // Take Deriv of iobj1 with respect to theta
void ddph(int iobj2,int iobj1); // Take Deriv of iobj1 with respect to phi
int irow2;
int ic2;
void gradient(int ior2,int iot2,int iop2,int iobj);
void curl(int ior2,int iot2,int iop2,int ior1,int iot1,int iop1);
void divergence(int iobj2,int ior1,int iot1,int iop1);
void cleanup(int iobj);
void multiply(int iobj3,int iobj2,int iobj1);
void add(int iobj3,int iobj2,int iobj1);
void subtract(int iobj3,int iobj2,int iobj1);
void legendre(int iobj1,int ordern);
void sphbessel(int iobj1,int orderL,int imater,int ipolar,int wavetype);
void copyobj(int iobj2,int iobj1);
void stress(int sigrr,int sigrth,int sigrph,int ur,int uth,int uph,int imater);
void multbyr(int iobj);
void divbyr(int iobj);
void divbyrsintheta(int iobj);
void multby(int iobj,double r1);
void multbyc(int iobj,complex cr1);
void settozero(int iobj);
void multbyA(int iobj);
void multbyB(int iobj);
void multbyC(int iobj);
void multbyD(int iobj);
void multbyE(int iobj);
void multbyF(int iobj);
int ngetm(int i1,int i2);
void nsetm(int i1,int i2,int i3);

double rho[NMATER],rho3=0.00001;
double cl3=3000;
double ct3=1000;
double css[2][NMATER];
double beta,omega,partrad_nm;
double klp,klm,kl3;
double ktp,ktm,kt3;
void evaluateobj(int iobj2,int iobj1,double partrad_nm);
complex factorofD(int iobj1);
complex factorofC(int iobj1);
complex factorofB(int iobj1);
complex factorofA(int iobj1);
complex factorofconst(int iobj1);
int mosttermseveradd=0;
int mosttermseverddr=0;
int mosttermsevercle=0;
void par(char *ch1);
void mat(char *ch1);
void makescr(void);

void main(void)
{ // main

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

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

double sfx=60,sfx2=10.5;
int bx0=24,by0=480-40;
int i77;

setcolor(9);
line(bx0,0,bx0,by0);
for (i77=5;i77<(bx0+320)/sfx2;i77=i77+5)
{ //
setcolor(8);
line(bx0+sfx2*i77,0,bx0+sfx2*i77,by0);
} //
setcolor(9);
line(bx0,by0,bx0+70*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]=1; 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[1]=3.6 ; css[0][1]=3.43; css[1][1]=2.14; mat("GeO2"); // GeO2
//rho[1]=2.2 ; css[0][1]=5.95; css[1][1]=3.76; mat("SiO2");
//rho[1]=2.33; css[0][1]=9.0; css[1][1]=5.2; mat("Si"); // rough avg.
rho[1]=0.5 ; css[0][1]=1; css[1][1]=0.3; mat("imag. mat'l.");
//rho[1]=1 ; css[0][1]=2; css[1][1]=0.7; mat("imag. mat'l.");
//rho[1]=1.5 ; css[0][1]=3; css[1][1]=1; mat("imag. 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]=0.25; css[0][1]=0.25*3.00; css[1][1]=0.25*1.732; 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.

rho[0]=rho[0]*1000;
rho[1]=rho[1]*1000;
css[0][0]=css[0][0]*1000;
css[1][0]=css[1][0]*1000;
css[0][1]=css[0][1]*1000;
css[1][1]=css[1][1]*1000;
gotoxy(5,15);cout<<rho[0]*0.001<<" g/cc cl="<<css[0][0]<<" ct="<<css[1][0]<<" m/s";
gotoxy(5,17);cout<<rho[1]*0.001<<" g/cc cl="<<css[0][1]<<" ct="<<css[1][1]<<" m/s";


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,1);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,2);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<<" ";

// Basic object is a series of terms.
int iobj;
for (iobj=0;iobj<NOBJ;iobj++)
{ // ..
bt[iobj]=-1; nt[iobj]=-1; mnt[iobj]=0;
// nt[]=-1 signifies variable not in use.
} // ..

int it1;
int ifac;
int io;
bt[0]=0;mnt[0]=140; if (bt[0]+mnt[0]>=NTRM) {cout<<"err 1110";exit(0);}
for (io=1;io<NOBJ;io++)
{ //
bt[io]=bt[io-1]+mnt[io-1];
mnt[io]=140;
if (bt[io]+mnt[io]>=NTRM) {cout<<"err 111 >NTRM "<<io;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

orderL=0;
gotoxy(5+22,18);cout<<"Spheroidal";
gotoxy(5+21,19);cout<<"Order l = "<<orderL<<" ";

/*
sphbessel(2,orderL,PARTIC,LONGIT,First_kind);
legendre(1,orderL);
multiply(0,1,2);
nt[1]=-1; // release variables
nt[2]=-1;
multbyA(0);
gradient(1,2,3,0);
nt[0]=-1;
*/
/*
ic2=1; irow2=4;
iobj=1;
gotoxy(ic2,irow2-1);cout<<" ";
gotoxy(ic2,irow2);  cout<<" ";
for (it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.. it1
dispterm(it1,ic2,irow2); irow2=irow2+2;
} //..
if (nt[iobj]==0) {gotoxy(ic2,irow2);cout<<"0";irow2=irow2+2;}
if (nt[iobj]<0)  {gotoxy(ic2,irow2);cout<<"?????";irow2=irow2+2;}
*/
/*
ic2=40; irow2=12;
iobj=2;
gotoxy(ic2,irow2-1);cout<<" ";
gotoxy(ic2,irow2);  cout<<" ";
for (it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.. it1
dispterm(it1,ic2,irow2); irow2=irow2+2;
} //..
if (nt[iobj]==0) {gotoxy(ic2,irow2);cout<<"0";irow2=irow2+2;}
if (nt[iobj]<0)  {gotoxy(ic2,irow2);cout<<"?????";irow2=irow2+2;}
*/

/*
sphbessel(4,orderL,PARTIC,TRANSV,First_kind);
legendre(5,orderL);
multiply(6,4,5);
nt[4]=-1; // release
nt[5]=-1;
multbyr(6);
multbyA(6);
settozero(7);
settozero(8);
curl(1,2,3,6,7,8);
nt[6]=-1; // release
nt[7]=-1;
nt[8]=-1;
*/
//curl(6,7,8,9,10,11);
//nt[9]=-1;
//nt[10]=-1;
//nt[11]=-1;

sphbessel(0,orderL,PARTIC,LONGIT,First_kind);
multbyA(0);
gradient(1,2,3,0);
nt[0]=-1;

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

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

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

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

sphbessel(8,orderL,MATRIX,LONGIT,Out_going);
multbyB(8);
sphbessel(11,orderL,MATRIX,LONGIT,In_going);
add(8,8,11);
nt[11]=-1;
//legendre(9,orderL);
//multiply(10,8,9);
//nt[8]=-1;
//nt[9]=-1;
//multbyr(10);
//settozero(8);
//settozero(9);
//curl(5,6,7,10,8,9);
//nt[8]=-1;
//nt[9]=-1;
//nt[10]=-1;
gradient(5,6,7,8);
nt[8]=-1;

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

getch();
return;
*/

//gotoxy(1,25);cout<<"s ";

/*
add(5,5,8);
add(6,6,9);
add(7,7,10);
nt[8]=-1;
nt[9]=-1;
nt[10]=-1;
*/

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

stress(8,9,10,1,2,3,PARTIC);
stress(11,12,13,5,6,7,MATRIX);
//nt[0]=-1;
//nt[1]=-1;

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


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

getch();
return;
*/
//gotoxy(1,25);cout<<"u ";

subtract(0,1,5); // r continuity
nt[1]=-1;
nt[5]=-1;
/*
subtract(1,2,6); // theta continuity
nt[2]=-1;
nt[6]=-1;
*/
/*
subtract(2,8,11); // rr force balance
nt[8]=-1;
nt[11]=-1;
*/
subtract(1,8,11); // r-r force balance
nt[8]=-1;
nt[11]=-1;

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

delay(400);

beta=0.1;

/*
gotoxy(1,1);cout<<"Number of terms for objects:";
for (i77=0;i77<24;i77++)
{
gotoxy(1,i77+2);cout<<i77<<". "<<nt[i77]<<" ";
}
gotoxy(40,10);cout<<"add: "<<mosttermseveradd<<" terms";
gotoxy(40,11);cout<<"ddr: "<<mosttermseverddr<<" terms";
gotoxy(40,12);cout<<"cle: "<<mosttermsevercle<<" terms";
*/

setcolor(11);
settextstyle(2,0,6);
outtextxy(bx0+0*sfx2-3, 442,"0");
outtextxy(bx0+5*sfx2-9,442,"5");
outtextxy(bx0+10*sfx2-9,442,"10");
outtextxy(bx0+15*sfx2-9,442,"15");
outtextxy(bx0+20*sfx2-9,442,"20");
outtextxy(bx0+25*sfx2-9,442,"25");
outtextxy(bx0+30*sfx2-9,442,"30");
outtextxy(bx0+40*sfx2-9,442,"40");
outtextxy(bx0+50*sfx2-9,442,"50");
outtextxy(bx0+60*sfx2-9,442,"60");
outtextxy(bx0+70*sfx2-9,442,"70");
settextstyle(2,0,6);
outtextxy(bx0+50+2,461,"Raman shift (cm  )");
settextstyle(2,0,5);
outtextxy(bx0+195+1,461-5,"-1");
/*
settextstyle(2,0,6);
outtextxy(bx0+0*sfx2-3,460,"0");
outtextxy(bx0+10*sfx2-9,460,"10");
outtextxy(bx0+20*sfx2-9,460,"20");
outtextxy(bx0+30*sfx2-9,460,"30");
outtextxy(bx0+35*sfx2-9,460,"cm");
settextstyle(2,0,5);
outtextxy(8+bx0+35*sfx2,460-6,"-1");
*/
settextstyle(2,1,6);
setcolor(13);
outtextxy(-1,190,"abs(A) (long)");
//setcolor(11);
//outtextxy(-1,330,"abs(B) (tran)");

partrad_nm=0.5*11.8; // Particle radius (arbitrary)
gotoxy(5,13);cout<<partrad_nm<<" nm radius";
gotoxy(26,13);cout<<2*partrad_nm<<" nm diam";

LB555:
if (beta<100.0*24.5*(0.30*2.0*3.14159*partrad_nm)/css[1][0])
{ //..
beta=beta+0.01;
} //..
else
{ //
makescr();
return;
} //
//gotoxy(1,1);cout<<"beta "<<beta<<" ";
omega=beta*css[1][0]/partrad_nm;
klp=omega/css[0][0];
ktp=omega/css[1][0];
if (NMATER>=2)
{ //
klm=omega/css[0][1];
ktm=omega/css[1][1];
} //
if (NMATER>=3)
{ //
kl3=omega/cl3;
kt3=omega/ct3;
} //

evaluateobj(4,0,partrad_nm);
evaluateobj(5,1,partrad_nm);

cleanup(4);
cleanup(5);

complex crb4,c3,c4,c5;
/*
crb4=factorofD(4);
c5=factorofD(5);
if (crb4==0) {cout<<"* crb4==0 1 *";getch();exit(0);}
crb4=-c5/crb4;
multbyc(4,crb4);
add(5,5,4);

crb4=factorofD(4);
c5=factorofD(6);
if (crb4==0) {cout<<"* crb4==0 2 *";getch();exit(0);}
crb4=-c5/crb4;
multbyc(4,crb4);
add(6,6,4);

crb4=factorofD(4);
c5=factorofD(7);
if (crb4==0) {cout<<"* crb4==0 3 *";getch();exit(0);}
crb4=-c5/crb4;
multbyc(4,crb4);
add(7,7,4);

crb4=factorofC(5);
c5=factorofC(6);
if (crb4==0) {cout<<"* crb4==0 4 *";getch();exit(0);}
crb4=-c5/crb4;
multbyc(5,crb4);
add(6,6,5);

crb4=factorofC(5);
c5=factorofC(7);
if (crb4==0) {cout<<"* crb4==0 5 *";getch();exit(0);}
crb4=-c5/crb4;
multbyc(5,crb4);
add(7,7,5);
*/

crb4=factorofB(4);
c5=factorofB(5);
if (crb4==0) {cout<<"* crb4==0 6 *";getch();exit(0);}
crb4=-c5/crb4;
multbyc(4,crb4);
add(5,5,4);
complex c51,c52,c53,c31,c3B,c3A;
c51=factorofA(4);
c52=factorofB(4);
c53=factorofconst(4);
// A*c51 + B*c52 + c53 = 0;
// B = -(A*c51 + c53)/c52 ;

c31=factorofA(5);
c3=-1.0/c31;
multbyc(5,c3);
double r67A,r67b,r67br,r67bi;
c3A=factorofconst(5);

c3B = -(c3A*c51+c53)/c52;

r67A=abs(c3A);
//r67A=r67A*beta; // re-scaling
r67b= fabs(abs(c3B));
r67br=fabs(real(c3B));
r67bi=fabs(imag(c3B));
// There needs to be appropriate frequency normalization, based
// on the density of states of phonons, assuming thermal
// occupation of all states.

float wavenum;
wavenum = 0.01*beta*css[1][0]/(0.3*2.0*3.14159*partrad_nm);

double sfy=60;

putpixel(bx0+sfx2*wavenum,by0-r67A*sfy,13);
setcolor(13);
  circle(bx0+sfx2*wavenum,by0-r67A*sfy,1);

//putpixel(bx0+sfx*beta,479-(1.0)*sfy,8);
//putpixel(bx0+sfx*beta,479-r67A*sfy*0.2,13);
//putpixel(bx0+sfx*beta,479-r67b*sfy,12-8);
//putpixel(bx0+sfx*beta,479-r67br*sfy,11-8);
//putpixel(bx0+sfx*beta,479-r67bi*sfy,10-8);

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

getch();
return;

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

} // main

void dispterm(int it2,int ic,int ir)
{ // dispterm
if (it2<0)     {cout<<"Error 10 ";exit(0);}
if (it2>=NTRM) {cout<<"Error 11 ";exit(0);}
if (ic<1)      {cout<<"Error 21 ";exit(0);}
if (ic>70)     {cout<<"Error 22 ";exit(0);}
if (ir<2)      {cout<<"Error 31 - 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<NFAC;ifac++)
{ //
if (ngetm(ifac,it2)!=0) anyfactors=1;
} //
if ((cr[it2]==1)&&(ci[it2]==0)&&(anyfactors==0))
{ //
cout<<"1";ic=ic+1;
} //

int ivar2;

for (ivar2=0;ivar2<2*NMATER;ivar2++)
{ // ivar2
if (ngetm(5+4*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<<"error 3553 ";getch();exit(0);}
ic=ic+3;
if (ngetm(5+4*ivar2,it2)!=1)
{ //..
gotoxy(ic,ir-1);cout<<ngetm(5+4*ivar2,it2)<<"  ";ic=ic+1;
} //..
gotoxy(ic,ir);cout<<" ";ic=ic+1;
} // kLp
} // ivar2

for (ivar2=-1;ivar2<2*NMATER;ivar2++)
{ // ivar2
if (ngetm(6+4*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<<"error 3554 ";getch();exit(0);}
gotoxy(ic,ir-1);if (ngetm(6+4*ivar2,it2)!=1) {cout<<ngetm(6+4*ivar2,it2)<<" ";ic=ic+1;}
gotoxy(ic,ir);cout<<" ";ic=ic+1;
} //..
} // ivar2

for (ivar2=-1;ivar2<2*NMATER;ivar2++)
{ // ivar2
if (ngetm(8+4*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<<"error 3555 ";getch();exit(0);}
gotoxy(ic,ir-1);if (ngetm(8+4*ivar2,it2)!=1) {cout<<ngetm(8+4*ivar2,it2)<<" ";ic=ic+1;}
gotoxy(ic,ir);cout<<" ";ic=ic+1;
} //..
} // ivar2

if (ngetm(3,it2)!=0)
{ //
gotoxy(ic,ir);cout<<"r  ";ic=ic+1;
if (ngetm(3,it2)!=1)
{ //..
gotoxy(ic,ir-1);cout<<ngetm(3,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(7+4*ivar,it2)!=0)
{ // A
gotoxy(ic,ir);cout<<sch1;ic=ic+1;
if (ngetm(7+4*ivar,it2)!=1)
{ //..
gotoxy(ic,ir-1);cout<<ngetm(7+4*ivar,it2)<<" ";ic=ic+1;
} //..
ic=ic+1;
} // A
} // ivar

return;
} // dispterm

void ddr(int iobj2,int iobj1)
{ // ddr
if (iobj2<0)       {cout<<"* Error ddr 1 ";getch();exit(0);}
if (iobj2>=NOBJ)   {cout<<"* Error ddr 2 ";getch();exit(0);}
if (iobj1<0)       {cout<<"* Error ddr 3 ";getch();exit(0);}
if (iobj1>=NOBJ)   {cout<<"* Error ddr 4 ";getch();exit(0);}
if (mnt[iobj1]<=0) {cout<<"Error ddr 5 mnt<0 ";getch();exit(0);}
if (mnt[iobj2]<=0) {cout<<"Error ddr 6 mnt<0";getch();exit(0);}
if (nt[iobj1]<0)   {cout<<"* Error ddr 7 nt<0";getch();exit(0);}
// Take derivative with respect to r:
int it2;
int it1;
it2=bt[iobj2];
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //.. it1

int i234;
for (i234=0;i234<=8*NMATER;i234=i234+4)
{ // i234

// Take derivative with respect to sin(r) term
if (it2>=bt[iobj2]+mnt[iobj2])
{gotoxy(1,23);cout<<"ddr error 8 >mnt Press a key";
 gotoxy(1,24);cout<<"obj="<<iobj2<<" "<<iobj1<<" ";
 gotoxy(1,25);cout<<"nt="<<nt[iobj2]<<" "<<nt[iobj1]<<" it2-bt="<<it2-bt[iobj2]<<"   ";
getch();exit(0);}
if (it2>=NTRM) {cout<<"ddr error 9 >NTRM Press a key";getch();exit(0);}
cr[it2]=cr[it1]*ngetm(i234+2,it1);
ci[it2]=ci[it1]*ngetm(i234+2,it1);
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
if (i234>0) nsetm(i234+1,it2,ngetm(i234+1,it1)+1);
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++;

// derivative with respect to cos(r) term
if (it2>=bt[iobj2]+mnt[iobj2]) {cout<<"ddr error 10 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"ddr error 11 >NTRM Press a key";getch();exit(0);}
cr[it2]=-cr[it1]*ngetm(i234+4, it1);
ci[it2]=-ci[it1]*ngetm(i234+4, it1);
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac,it1));
} //
if (i234>0) nsetm(i234+1, it2, ngetm(i234+1, it1)+1);
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

int ifac;
// derivative with respect to r^n3 term
if (it2>=bt[iobj2]+mnt[iobj2]) {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(3, it1);
ci[it2]=ci[it1]*ngetm(3, it1);
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac, it2, ngetm(ifac,it1));
} //
nsetm(3, it2, ngetm(3, it1)-1); // r
if ((cr[it2]!=0)||(ci[it2]!=0)) it2++;
} //.. it1
nt[iobj2]=it2-bt[iobj2];
if (nt[iobj2]>mnt[iobj2]) {cout<<"* Err ddr nt>mnt";getch();exit(0);}
if (nt[iobj2]>mosttermseverddr) {mosttermseverddr=nt[iobj2];}
return;
} // ddr

void ddth(int iobj2,int iobj1)
{ // ddth
if (iobj1<0)       {cout<<"* Error ddth 1 iobj1";getch();exit(0);}
if (iobj1>=NOBJ)   {cout<<"* Error ddth 2 iobj1";getch();exit(0);}
if (iobj2<0)       {cout<<"* Error ddth 3 iobj2";getch();exit(0);}
if (iobj2>=NOBJ)   {cout<<"* Error ddth 4 iobj2";getch();exit(0);}
if (mnt[iobj1]<=0) {cout<<"* Error ddth 5 mnt<0";getch();exit(0);}
if (mnt[iobj2]<=0) {cout<<"* Error ddth 6 mnt<0";getch();exit(0);}
if (nt[iobj1]<0)   {cout<<"* Error ddth 7 nt< 0";getch();exit(0);}
int it2;
int it1;
it2=bt[iobj2];
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //.. it1
// Take derivative with respect to sin(th) term
if (it2>=bt[iobj2]+mnt[iobj2]) {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);
int ifac;
for (ifac=0;ifac<NFAC;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[iobj2]+mnt[iobj2]) {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<NFAC;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[iobj2]=it2-bt[iobj2];
if (nt[iobj2]>mnt[iobj2]) {cout<<"* Err ddth nt>mnt";getch();exit(0);}
return;
} // ddth

void ddph(int iobj2,int iobj1)
{ // ddph
if (iobj2<0)       {cout<<"* Error ddph 1 ";getch();exit(0);}
if (iobj2>=NOBJ)   {cout<<"* Error ddph 2 ";getch();exit(0);}
if (iobj1<0)       {cout<<"* Error ddph 3 ";getch();exit(0);}
if (iobj1>=NOBJ)   {cout<<"* Error ddph 4 ";getch();exit(0);}
if (mnt[iobj1]<=0) {cout<<"Error ddph 5 mnt<0 ";getch();exit(0);}
if (mnt[iobj2]<=0) {cout<<"Error ddph 6 mnt<0";getch();exit(0);}
if (nt[iobj1]<0)   {cout<<"* Error ddph 7 nt<0";getch();exit(0);}
int it2;
int it1;
it2=bt[iobj2];
nt[iobj2]=it2-bt[iobj2];
return;
} // ddph

void gradient(int ior2,int iot2,int iop2,int iobj)
{ // gradient
if (nt[ior2 ]!=-1)    {cout<<"* Err gradie 0 var "<<ior2 <<" in use. ";getch();exit(0);}
if (nt[iot2 ]!=-1)    {cout<<"* Err gradie 0 var "<<iot2 <<" in use. ";getch();exit(0);}
if (nt[iop2 ]!=-1)    {cout<<"* Err gradie 0 var "<<iop2 <<" in use. ";getch();exit(0);}
int it1;
ddr(ior2,iobj);
ddth(iot2,iobj);
// Divide object by r
for (it1=bt[iot2];it1<bt[iot2]+nt[iot2];it1++)
{ //.. it1
nsetm(3, it1, ngetm(3, it1)-1);
} //.. it1
ddph(iop2,iobj);
// Divide object by r sin(theta)
for (it1=bt[iop2];it1<bt[iop2]+nt[iop2];it1++)
{ //.. it1
nsetm(3, it1, ngetm(3, it1)-1);
nsetm(0, it1, ngetm(0, it1)-1);
} //.. it1
return;
} // gradient

void divergence(int iobj2,int ior1,int iot1,int iop1)
{ // divergence
int it1,ioaux;
ioaux=NOBJ-1; // auxilliary object (temp variable)
if (nt[ioaux]!=-1) {cout<<"error div";getch();exit(0);}

// mult by r^2
for (it1=bt[ior1];it1<bt[ior1]+nt[ior1];it1++)
{ //.. it1
nsetm(3, it1, ngetm(3, it1)+2); // r
} //.. it1

ddr(ioaux,ior1);

// divide by r^2 (to restore)
for (it1=bt[ior1];it1<bt[ior1]+nt[ior1];it1++)
{ //.. it1
nsetm(3, it1, ngetm(3, it1)-2); // r
} //.. it1

// divide by r^2
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //.. it1
nsetm(3,it1,ngetm(3, it1)-2); // r
} //.. it1

int it2;
it2=bt[iobj2];
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //.. it1
if ((cr[it1]!=0)||(ci[it1]!=0))
{ //..
if (it2>=bt[iobj2]+mnt[iobj2]) {cout<<"ghj error 2 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"ghj error 4 >NTRM Press a key";getch();exit(0);}
cr[it2]=cr[it1]; ci[it2]=ci[it1];
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac, it2, ngetm(ifac, it1));
} //
it2=it2+1;
} //
} //.. it1
nt[iobj2]=it2-bt[iobj2];

// mult by sin(theta)
for (it1=bt[iot1];it1<bt[iot1]+nt[iot1];it1++)
{ //.. it1
nsetm(0, it1, ngetm(0, it1)+1); // sin(theta)
} //.. it1

ddth(ioaux,iot1);

// divide by sin(theta) to restore
for (it1=bt[iot1];it1<bt[iot1]+nt[iot1];it1++)
{ //.. it1
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
} //.. it1

// divide by r sin(theta)
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //.. it1
nsetm(3, it1, ngetm(3, it1)-1); // r
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
} //.. it1

// Concatenate (add on)
it2=bt[iobj2]+nt[iobj2];
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //.. it1
if ((cr[it1]!=0)||(ci[it1]!=0))
{ //..
if (it2>=bt[iobj2]+mnt[iobj2]) {cout<<"sdf error 2 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"dfg error 4 >NTRM Press a key";getch();exit(0);}
cr[it2]=cr[it1];
ci[it2]=ci[it1];
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac,it2,ngetm(ifac, it1));
} //
it2=it2+1;
} //..
} //.. it1
nt[iobj2]=it2-bt[iobj2];

ddph(ioaux,iop1);

// divide by r sin theta
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //.. it1
nsetm(3, it1, ngetm(3, it1)-1); // r
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
} //.. it1

// Concatenate (add on)
it2=bt[iobj2]+nt[iobj2];
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //.. it1
if ((cr[it1]!=0)||(ci[it1]!=0))
{ //..
if (it2>=bt[iobj2]+mnt[iobj2]) {cout<<"asd error 2 >mnt Press a key";getch();exit(0);}
if (it2>=NTRM) {cout<<"asd error 4 >NTRM Press a key";getch();exit(0);}
cr[it2]=cr[it1];
ci[it2]=ci[it1];
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac, it2, ngetm(ifac, it1));
} //
it2=it2+1;
} //..
} //.. it1
nt[iobj2]=it2-bt[iobj2];
if (nt[iobj2]>mnt[iobj2]) {cout<<"* Err diverg nt>mnt";getch();exit(0);}

nt[ioaux]=-1;

return;
} // divergence

void cleanup(int iobj)
{ // cleanup
if (nt[iobj]>mosttermsevercle) {mosttermsevercle=nt[iobj];}
int it1,it2,it3;
int flag1;

int i234;
for (i234=0;i234<=2*NMATER;i234=i234+4)
{ // i234
//  Look for reducible factors of cos(r)^2, cos(kLp r), cos(kTp r):
flag1=0;
doagain:
it3=bt[iobj]+nt[iobj];
for (it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.. it1
if (ngetm(i234+4, it1)>=2)
{ //.
if (it3>=bt[iobj]+mnt[iobj]) {cout<<"iii error 2 >mnt Press a key";getch();exit(0);}
if (it3>=NTRM)               {cout<<"kkk error 4 >NTRM Press a key";getch();exit(0);}
flag1=1;
nsetm(i234+4, it1, ngetm(i234+4, it1)-2); // cos(r)
it2=bt[iobj]+nt[iobj];nt[iobj]=nt[iobj]+1;
cr[it2]=-1*cr[it1];
ci[it2]=-1*ci[it1];
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac, it2, ngetm(ifac, it1));
} //
nsetm(i234+2, it2, ngetm(i234+2, it1)+2); // sin(r)
} //.
} // it1
if (flag1==1) goto doagain;
} // i234

// Next, look for reducible factors of cos(theta)^2
flag1=0;
doagain2:
it3=bt[iobj]+nt[iobj];
for (it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.. it1
if (ngetm(1, it1)>=2)
{ //.
if (it3>=bt[iobj]+mnt[iobj]) {cout<<"ggg error 2 >mnt Press a key";getch();exit(0);}
if (it3>=NTRM) {cout<<"lll error 4 >NTRM Press a key";getch();exit(0);}
nsetm(1, it1, ngetm(1, it1)-2); // divide by cos(theta)^2
it2=bt[iobj]+nt[iobj]; nt[iobj]=nt[iobj]+1;
cr[it2]=-1*cr[it1]; ci[it2]=-1*ci[it1];
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac, it2, ngetm(ifac,it1));
} //
nsetm(0, it2, ngetm(0,it1)+2); // sin(theta)
} //.
} //.. it1
if (flag1==1) goto doagain2;

// Now add up terms of equal powers in all factors:
//  (thus eliminating extra terms)
startagain:
for (it1=bt[iobj];it1<bt[iobj]+nt[iobj]-1;it1++)
{ //.. it1
for (it2=it1+1;it2<bt[iobj]+nt[iobj];it2++)
{ //... it2
int flagsame;
flagsame=1;
int ifac;
for (ifac=0;ifac<NFAC;ifac++)
{ //
if (ngetm(ifac,it1)!=ngetm(ifac,it2)) flagsame=0;
} //
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[iobj]+nt[iobj]-1; // last term in object
if (it2!=it3)
{ //.. not last term
cr[it2]=cr[it3]; ci[it2]=ci[it3];
for (ifac=0;ifac<NFAC;ifac++)
{ //.
nsetm(ifac, it2, ngetm(ifac, it3));
} //.
} //.. not last term
nt[iobj]=nt[iobj]-1;
it3=bt[iobj]+nt[iobj]-1; // last term in object
if ((cr[it1]==0)&&(ci[it1]==0))
{ //..
if (it1!=it3)
{ //.. not last term
cr[it1]=cr[it3]; ci[it1]=ci[it3];
for (ifac=0;ifac<NFAC;ifac++)
{ //.
nsetm(ifac, it1, ngetm(ifac, it3));
} //.
} //.. not last term
nt[iobj]=nt[iobj]-1;
} //..
goto startagain;
} //. found two terms
} //... it2
} //.. it1

return;
} // cleanup

void multiply(int iobj3,int iobj2,int iobj1)
{ // multiply
if (nt[iobj3]!=-1)    {cout<<"* Err multip 0 var "<<iobj3<<" in use. ";getch();exit(0);}
if ((iobj3==iobj1)||(iobj3==iobj2))  {cout<<"multiply error 1  Press a key";getch();exit(0);}
if (iobj3<0)     {cout<<"* Error mult 2 ";getch();exit(0);}
if (iobj3>=NOBJ) {cout<<"* Error mult 3 ";getch();exit(0);}
if (iobj2<0)     {cout<<"* Error mult 4 ";getch();exit(0);}
if (iobj2>=NOBJ) {cout<<"* Error mult 5 ";getch();exit(0);}
if (iobj1<0)     {cout<<"* Error mult 6 ";getch();exit(0);}
if (iobj1>=NOBJ) {cout<<"* Error mult 7 ";getch();exit(0);}
if (mnt[iobj1]<=0) {cout<<"Error mult 8 mnt<0 ";getch();exit(0);}
if (mnt[iobj2]<=0) {cout<<"Error mult 9 mnt<0";getch();exit(0);}
if (mnt[iobj3]<=0) {cout<<"Error mult 10 mnt<0";getch();exit(0);}
int it1,it2,it3,ifac;
it3=bt[iobj3];
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ // it1
for (it2=bt[iobj2];it2<bt[iobj2]+nt[iobj2];it2++)
{ // it2
if (it3>=bt[iobj3]+mnt[iobj3]) {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<NFAC;ifac++)
{ //.
nsetm(ifac, it3, ngetm(ifac,it1)+ngetm(ifac, it2));
} //.
it3++;
} // it2
} // it1
nt[iobj3]=it3-bt[iobj3];
if (nt[iobj3]>mnt[iobj3]) {cout<<"Error mult nt>mnt";getch();exit(0);}
cleanup(iobj3);
return;
} // multiply

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

void copyobj(int iobj2,int iobj1)
{ // copyobj
if (iobj2<0)     {cout<<"* Error copyobj 4 ";getch();exit(0);}
if (iobj2>=NOBJ) {cout<<"* Error copyobj 5 ";getch();exit(0);}
if (iobj1<0)     {cout<<"* Error copyobj 6 ";getch();exit(0);}
if (iobj1>=NOBJ) {cout<<"* Error copyobj 7 ";getch();exit(0);}
if (mnt[iobj1]<=0) {cout<<"Error copyobj 8 mnt<0 ";getch();exit(0);}
if (mnt[iobj2]<=0)
{
gotoxy(1,22);cout<<"Error copyobj 9 mnt["<<iobj2<<"]="<<mnt[iobj2]<<" ";
getch();exit(0);
}
if (mnt[iobj2]<nt[iobj1]) {cout<<"Error copyobj 10 mnt<nt";getch();exit(0);}
int it1,ifac,it2;
it2=bt[iobj2];
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //..
cr[it2]=cr[it1]; ci[it2]=ci[it1];
for (ifac=0;ifac<NFAC;ifac++) nsetm(ifac, it2, ngetm(ifac, it1));
it2++;
} //..
nt[iobj2]=nt[iobj1];
if (nt[iobj2]>mnt[iobj2])
{
gotoxy(1,22);cout<<"Error copyobj 11 nt>mnt: "<<iobj2<<" "<<nt[iobj2]<<" > "<<mnt[iobj2]<<" ";
getch();exit(0);
}
return;
} // copyobj

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   0 var "<<iot2 <<" in use. ";getch();exit(0);}
if (nt[iop2 ]!=-1)    {cout<<"* Err curl   0 var "<<iop2 <<" in use. ";getch();exit(0);}
// This assumes that there is no phi dependence of any
// of the components.  One aux. object is needed.
int ioaux1; ioaux1=NOBJ-1;
if (nt[ioaux1]!=-1) {cout<<"curl error 1 var in use";getch();exit(0);}
int it1;
for (it1=bt[iop1];it1<bt[iop1]+nt[iop1];it1++)
{ //
nsetm(0, it1, ngetm(0, it1)+1); // sin(theta)
} //
ddth(ior2,iop1);
// Restore iop1:
for (it1=bt[iop1];it1<bt[iop1]+nt[iop1];it1++)
{ //
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
} //
for (it1=bt[ior2];it1<bt[ior2]+nt[ior2];it1++)
{ //
nsetm(3, it1, ngetm(3, it1)-1); // r
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
} //
//------------------------------
for (it1=bt[iop1];it1<bt[iop1]+nt[iop1];it1++)
{ //
nsetm(3, it1, ngetm(3, it1)+1); // r
} //
ddr(iot2,iop1);
// Restore iop1:
for (it1=bt[iop1];it1<bt[iop1]+nt[iop1];it1++)
{ //
nsetm(3, it1, ngetm(3,it1)-1); // r
} //
for (it1=bt[iot2];it1<bt[iot2]+nt[iot2];it1++)
{ //.
cr[it1]=cr[it1]*(-1);
ci[it1]=ci[it1]*(-1);
nsetm(3, it1, ngetm(3,it1)-1); // r
} //.
//------------------------------
for (it1=bt[iot1];it1<bt[iot1]+nt[iot1];it1++)
{ //
nsetm(3, it1, ngetm(3, it1)+1); // r
} //
ddr(iop2,iot1);
for (it1=bt[iot1];it1<bt[iot1]+nt[iot1];it1++)
{ //
nsetm(3, it1, ngetm(3, it1)-1); // r
} //
ddth(ioaux1,ior1);
for (it1=bt[ioaux1];it1<bt[ioaux1]+nt[ioaux1];it1++)
{ //.
cr[it1]=cr[it1]*(-1);
ci[it1]=ci[it1]*(-1);
} //.
add(iop2,iop2,ioaux1);
for (it1=bt[iop2];it1<bt[iop2]+nt[iop2];it1++)
{ //.
nsetm(3, it1, ngetm(3, it1)-1); // r
} //.
nt[ioaux1]=-1; // release aux var
return;
} // curl

void add(int iobj3,int iobj2,int iobj1)
{ // add
//if (nt[iobj3]!=-1)   {cout<<"* Err add 0 var "<<iobj3<<" in use. ";getch();exit(0);}
if (nt[iobj1]+nt[iobj2]>mnt[iobj3]) {cout<<"* Error add 0 ";getch();exit(0);}
if (iobj1<0)     {cout<<"* Error add 1 ";getch();exit(0);}
if (iobj1>=NOBJ) {cout<<"* Error add 2 ";getch();exit(0);}
if (iobj2<0)     {cout<<"* Error add 3 ";getch();exit(0);}
if (iobj2>=NOBJ) {cout<<"* Error add 4 ";getch();exit(0);}
if (iobj3<0)     {cout<<"* Error add 5 ";getch();exit(0);}
if (iobj3>=NOBJ) {cout<<"* Error add 6 ";getch();exit(0);}
// It is okay for any of these three objects to be the same.
int i3,it2,it3;
if (iobj2==iobj3) { i3=iobj2; iobj2=iobj1; iobj1=i3; }
copyobj(iobj3,iobj1);
it3=bt[iobj3]+nt[iobj3];
for (it2=bt[iobj2];it2<bt[iobj2]+nt[iobj2];it2++)
{ // it2
if (it3>=bt[iobj3]+mnt[iobj3])
{ //..
gotoxy(1,21);cout<<"* add error 7 >mnt Press a key";
gotoxy(1,22);cout<<it3<<" "<<bt[iobj3]<<" "<<mnt[iobj3]<<" ";
gotoxy(1,23);cout<<"* obj: "<<iobj3<<" "<<iobj2<<" "<<iobj1<<" *";
gotoxy(1,24);cout<<"* nt: "<<nt[iobj3]<<", "<<nt[iobj2]<<", "<<nt[iobj1]<<" terms *";
gotoxy(1,25);cout<<"* mnt: "<<mnt[iobj3]<<", "<<mnt[iobj2]<<", "<<mnt[iobj1]<<" terms *";
getch();exit(0);
} //..
if (it3>=NTRM) {cout<<"add error 8 >NTRM Press a key";getch();exit(0);}
int ifac;
cr[it3]=cr[it2]; ci[it3]=ci[it2];
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac,it3,ngetm(ifac,it2));
} //
it3=it3+1;
} // it2
nt[iobj3]=it3-bt[iobj3];
if (nt[iobj3]>mnt[iobj3]) {cout<<"93854";getch();exit(0);}
if (nt[iobj3]>mosttermseveradd) {mosttermseveradd=nt[iobj3];}
cleanup(iobj3);
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;
int okay1=0;
if (imater==PARTIC)
{ //..
shearmodmu=rho[0]*css[1][0]*css[1][0];
lambdalame=rho[0]*css[0][0]*css[0][0]-2*shearmodmu;
okay1=1;
} //..
if (imater==MATRIX)
{ //..
shearmodmu=rho[1]*css[1][1]*css[1][1];
lambdalame=rho[1]*css[0][1]*css[0][1]-2*shearmodmu;
okay1=1;
} //..
if (imater==MATER3)
{ //..
shearmodmu=rho3*ct3*ct3;
lambdalame=rho3*cl3*cl3-2*shearmodmu;
okay1=1;
} //..
if (okay1==0) {cout<<"stress err";getch();exit(0);}

// Reference http://www.rochester.edu/College/ME/gans/ME445/stress.pdf
// c:/cofrest/stress.pdf
int err,ethth,ephph,lamdilat;
int ioaux1;
ioaux1=NOBJ-1;
err=NOBJ-2;
ethth=NOBJ-3;
ephph=NOBJ-4;
lamdilat=NOBJ-5;
if (nt[ioaux1]!=-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 (ioaux1==sigrr) {cout<<"* err stress 17 var conflic.";getch();exit(0);}
if (ioaux1==sigrth) {cout<<"* err stress 18 var conflic.";getch();exit(0);}
if (ioaux1==sigrph) {cout<<"* err stress 19 var conflic.";getch();exit(0);}
int it1;
//-----------------
ddr(err,ur);
//-----------------
ddth(ethth,uth);
add(ethth,ethth,ur);
divbyr(ethth);
//-----------------
copyobj(ephph,uth);
for (it1=bt[ephph];it1<bt[ephph]+nt[ephph];it1++)
{ //.
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
nsetm(1, it1, ngetm(1, it1)+1); // cos(theta)
} //.
add(ephph,ephph,ur);
divbyr(ephph);
//-----------------------
ddr(sigrph,uph);
copyobj(ioaux1,uph);
multby(ioaux1,-1.0);
divbyr(ioaux1);
add(sigrph,sigrph,ioaux1);
multby(sigrph,0.5);
//---------------------
ddth(sigrth,ur);
divbyr(sigrth);
ddr(ioaux1,uth);
add(sigrth,sigrth,ioaux1);
copyobj(ioaux1,uth);
multby(ioaux1,-1.0);
divbyr(ioaux1);
add(sigrth,sigrth,ioaux1);
multby(sigrth,0.5);
//------------------------
add(lamdilat,err,ethth);
add(lamdilat,lamdilat,ephph);
multby(lamdilat,lambdalame);

copyobj(sigrr,err);
multby(sigrr,2.0*shearmodmu);
add(sigrr,sigrr,lamdilat);

multby(sigrph,2.0*shearmodmu);
multby(sigrth,2.0*shearmodmu);

nt[ioaux1]=-1;
nt[err]=-1;
nt[ethth]=-1;
nt[ephph]=-1;
nt[lamdilat]=-1;

return;
} // strain

void multbyr(int iobj)
{ // multbyr
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.
nsetm(3, it1, ngetm(3, it1)+1); // r
} //.
return;
} // multbyr

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

void divbyr(int iobj)
{ // divbyr
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.
nsetm(3, it1, ngetm(3,it1)-1); // r
} //.
return;
} // divbyr

void divbyrsintheta(int iobj)
{ // divbyrsintheta
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++)
{ //.
nsetm(3, it1, ngetm(3, it1)-1); // r
nsetm(0, it1, ngetm(0, it1)-1); // sin(theta)
} //.
return;
} // divbyrsintheta

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

void multbyA(int iobj)
{ // multbyA
if (NFAC<7+1) {cout<<"error A NFAC";getch();exit(0);}
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++) nsetm(7, it1, ngetm(7, it1)+1);
return;
} // multbyA

void multbyB(int iobj)
{ // multbyB
if (NFAC<11+1) {cout<<"error B NFAC";getch();exit(0);}
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++) nsetm(11, it1, ngetm(11, it1)+1);
return;
} // multbyB

void multbyC(int iobj)
{ // multbyC
if (NMATER<2)  {cout<<"error C NMAT";getch();exit(0);}
if (NFAC<15+1) {cout<<"error C NFAC";getch();exit(0);}
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++) nsetm(15, it1, ngetm(15, it1)+1);
return;
} // multbyC

void multbyD(int iobj)
{ // multbyD
if (NMATER<2)  {cout<<"error D NMAT";getch();exit(0);}
if (NFAC<19+1) {cout<<"error D NFAC";getch();exit(0);}
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++) nsetm(19, it1, ngetm(19, it1)+1);
return;
} // multbyD

void multbyE(int iobj)
{ // multbyE
if (NMATER<3)  {cout<<"error E NMAT";getch();exit(0);}
if (NFAC<23+1) {cout<<"error E NFAC";getch();exit(0);}
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++) nsetm(23, it1, ngetm(23, it1)+1);
return;
} // multbyE

void multbyF(int iobj)
{ // multbyF
if (NMATER<3)  {cout<<"error F NMAT";getch();exit(0);}
if (NFAC<27+1) {cout<<"error F NFAC";getch();exit(0);}
for (int it1=bt[iobj];it1<bt[iobj]+nt[iobj];it1++) nsetm(27, it1,ngetm(27, it1)+1);
return;
} // multbyF

void subtract(int iobj3,int iobj2,int iobj1)
{ // subtract
if (nt[iobj3]!=-1) {cout<<"* Err subtra 0 var "<<iobj3<<" in use. ";getch();exit(0);}
if (iobj1<0)     {cout<<"* Error subtract 1 ";getch();exit(0);}
if (iobj1>=NOBJ) {cout<<"* Error subtract 2 ";getch();exit(0);}
if (iobj2<0)     {cout<<"* Error subtract 3 ";getch();exit(0);}
if (iobj2>=NOBJ) {cout<<"* Error subtract 4 ";getch();exit(0);}
if (iobj3<0)     {cout<<"* Error subtract 5 ";getch();exit(0);}
if (iobj3>=NOBJ) {cout<<"* Error subtract 6 ";getch();exit(0);}
// It is okay for any of these three objects to be the same.
int i3,it2,it3,signflag=-1;
int ifac;
if (iobj2==iobj3) { i3=iobj2; iobj2=iobj1; iobj1=i3; signflag*=-1; }
it3=bt[iobj3];
for (it2=bt[iobj1];it2<bt[iobj1]+nt[iobj1];it2++)
{ // it2
if (it3>=bt[iobj3]+mnt[iobj3])
{ //
gotoxy(1,22);cout<<"subtract error 7 >mnt Press a key";
gotoxy(1,23);cout<<"obj: "<<iobj3<<" "<<iobj2<<" "<<iobj1<<" ";
getch();exit(0);
} //
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<NFAC;ifac++)
{ //
nsetm(ifac, it3, ngetm(ifac,it2));
} //
it3=it3+1;
} //
signflag=signflag*(-1);
it3=bt[iobj3]+nt[iobj1];
for (it2=bt[iobj2];it2<bt[iobj2]+nt[iobj2];it2++)
{ // it2
if (it3>=bt[iobj3]+mnt[iobj3])
{
gotoxy(1,22);cout<<"subtract error 9 >mnt Press a key";
gotoxy(1,23);cout<<"obj: "<<iobj3<<" "<<iobj2<<" "<<iobj1<<" ";
getch();exit(0);
}
if (it3>=NTRM)                 {cout<<"subtract error 10>NTRM Press a key";getch();exit(0);}
cr[it3]=cr[it2]*signflag;
ci[it3]=ci[it2]*signflag;
for (ifac=0;ifac<NFAC;ifac++)
{ //
nsetm(ifac, it3, ngetm(ifac,it2));
} //
it3=it3+1;
} // it2
nt[iobj3]=it3-bt[iobj3];
cleanup(iobj3);
return;
} // subtract

void sphbessel(int iobj1,int orderL,int imater,int ipolar,int wavetype)
{ // sphbessel
if (nt[iobj1]!=-1)    {cout<<"* Err sphbes 0 var "<<iobj1<<" 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 (iobj1<0)          {cout<<"* Err sphbes 6 ";getch();exit(0);}
if (iobj1>=NOBJ)      {cout<<"* Err sphbes 7 ";getch();exit(0);}
if (mnt[iobj1]<=0)    {cout<<"* Err sphbes 8 mnt<0 ";getch();exit(0);}
int it1,ifac,ioaux,it2,iL;
ioaux=NOBJ-1;
if (nt[ioaux]!=-1)    {cout<<"* Err sphbess 9";getch();exit(0);}

it1=bt[ioaux];
for (ifac=0;ifac<NFAC;ifac++) nsetm(ifac, it1, 0);
cr[it1]=1; ci[it1]=0;
nsetm(3, it1, -1); // r
if ((wavetype==0)||(wavetype==2)||(wavetype==3))
{ // 0 = 1st
nsetm(2+4*kfactor, it1, 1); // sin(r)
if (kfactor>0) nsetm(1+4*kfactor, it1, -1);
} // 0 = 1st
if (wavetype==1)
{ // 1 = 2nd
nsetm(4+4*kfactor, it1, 1); // cos(r)
if (kfactor>0) nsetm(1+4*kfactor, it1, -1);
} // 1 = 2nd
nt[ioaux]=1;
if ((wavetype==2)||(wavetype==3))
{ // 3 = outgoing 2 = ingoing
it1=bt[ioaux]+1;
for (ifac=0;ifac<NFAC;ifac++) nsetm(ifac, it1, 0);
cr[it1]=0;
if (wavetype==3) ci[it1]=1; // outgoing
if (wavetype==2) ci[it1]=-1; // ingoing
nsetm(3, it1, -1); // r
nsetm(4+4*kfactor, it1, 1); // cos(r)
if (kfactor>0) nsetm(1+4*kfactor, it1, -1); // k
nt[ioaux]=2;
} // 3 = outgoing  2 = ingoing

for (iL=0;iL<orderL;iL++)
{ //. iL
ddr(iobj1,ioaux);
it1=bt[iobj1];
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //..
nsetm(3, it1, ngetm(3,it1)-1); // r
if (kfactor>0) nsetm(1+4*kfactor,it1,ngetm(1+4*kfactor,it1)-2); // k factor
} //..
copyobj(ioaux,iobj1);
} //. iL
for (iL=0;iL<orderL;iL++)
{ //. iL
it1=bt[ioaux];
for (it1=bt[ioaux];it1<bt[ioaux]+nt[ioaux];it1++)
{ //..
cr[it1]=cr[it1]*(-1); ci[it1]=ci[it1]*(-1);
nsetm(3, it1, ngetm(3,it1)+1); // r
if (kfactor>0) nsetm(1+4*kfactor,it1,ngetm(1+4*kfactor,it1)+1); // k factor
} //..
} //. iL
copyobj(iobj1,ioaux);

nt[ioaux]=-1;

return;
} // sphbessel

void evaluateobj(int iobj2,int iobj1,double partrad_nm)
{ // evaluateobj

if (mnt[iobj2]<nt[iobj1])
{gotoxy(1,20);cout<<"* error evalobj 1 mnt<nt";
gotoxy(1,21);cout<<"* obj: "<<iobj2<<" "<<iobj1<<" ";
gotoxy(1,22);cout<<"* mnt[iobj2]="<<mnt[iobj2]<<" nt[iobj1]="<<nt[iobj1]<<"  ";
getch();exit(0);}

double r1,r2,r3;
int it1,it2,ifac;
it1=bt[iobj1];
nt[iobj2]=nt[iobj1];
for (it2=bt[iobj2];it2<bt[iobj2]+nt[iobj2];it2++)
{ //. it2
cr[it2]=cr[it1]; ci[it2]=ci[it1];
int flag3;
for (ifac=0;ifac<NFAC;ifac++)
{ //..
nsetm(ifac, it2, ngetm(ifac,it1));
flag3=1;
if (ifac==7)  flag3=0;
if (ifac==11) flag3=0;
if (ifac==15) flag3=0;
if (ifac==19) flag3=0;
if (flag3==1) nsetm(ifac, it2, 0);
} //..
it1++;
} //. it2

it2=bt[iobj2];
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //.
if (ngetm(3,it1)!=0)
{ // r
r1=partrad_nm; r2=pow(r1,ngetm(3,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} // r
if (ngetm(5,it1)!=0)
{ //
r1=klp; r2=pow(r1,ngetm(5,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(9,it1)!=0)
{ //
r1=ktp; r2=pow(r1,ngetm(9,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(13,it1)!=0)
{ //
r1=klm; r2=pow(r1,ngetm(13,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(17,it1)!=0)
{ //
r1=ktm; r2=pow(r1,ngetm(17,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(6,it1)!=0)
{ //
r1=sin(klp*partrad_nm); r2=pow(r1,ngetm(6,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(8,it1)!=0)
{ //
r1=cos(klp*partrad_nm); r2=pow(r1,ngetm(8,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(10,it1)!=0)
{ //
r1=sin(ktp*partrad_nm); r2=pow(r1,ngetm(10,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(12,it1)!=0)
{ //
r1=cos(ktp*partrad_nm); r2=pow(r1,ngetm(12,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(14,it1)!=0)
{ //
r1=sin(klm*partrad_nm); r2=pow(r1,ngetm(14,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(16,it1)!=0)
{ //
r1=cos(klm*partrad_nm); r2=pow(r1,ngetm(16,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(18,it1)!=0)
{ //
r1=sin(ktm*partrad_nm); r2=pow(r1,ngetm(18,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
if (ngetm(20,it1)!=0)
{ //
r1=cos(ktm*partrad_nm); r2=pow(r1,ngetm(20,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} //
it2++;
} //.
return;
} // evaluateobj

complex factorofD(int iobj1)
{ //
complex cr3;
cr3=0.0;
int it1;
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //
if (ngetm(19,it1)!=0)
{
cr3=complex(cr[it1],ci[it1]);
}
} //
return cr3;
} //

complex factorofC(int iobj1)
{ //
complex cr3;
cr3=complex(0.0,0.0);
int it1;
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //
if (ngetm(15, it1)!=0)
{
cr3=complex(cr[it1],ci[it1]);
}
} //
return cr3;
} //

complex factorofB(int iobj1)
{ //
complex cr3;
cr3=0.0;
int it1;
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //
if (ngetm(11, it1)!=0)
{
cr3=complex(cr[it1],ci[it1]);
}
} //
return cr3;
} //

complex factorofA(int iobj1)
{ //
complex cr3;
cr3=0.0;
int it1;
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //
if (ngetm(7,it1)!=0)
{
cr3=complex(cr[it1],ci[it1]);
}
} //
return cr3;
} //

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

complex factorofconst(int iobj1)
{ //
complex cr3;
cr3=0.0;
int it1;
for (it1=bt[iobj1];it1<bt[iobj1]+nt[iobj1];it1++)
{ //
if ((ngetm(11,it1)==0)&&(ngetm(7,it1)==0)&&(ngetm(15,it1)==0)&&(ngetm(19,it1)==0))
{
cr3=complex(cr[it1],ci[it1]);
}
} //
return cr3;
} //

int ngetm(int i1,int i2)
{ // ngetm
if (i1<0    ) {cout<<"err 1353";getch();exit(0);}
if (i2<0    ) {cout<<"err 1536";getch();exit(0);}
if (i1>=NFAC) {cout<<"err 5353";getch();exit(0);}
if (i2>=NTRM) {cout<<"err 5536";getch();exit(0);}
int i3;
if (i2<ntrm0) i3=nm0[i1][i2];
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)];
return i3;
} // ngetm

void nsetm(int i1,int i2,int i3)
{ // nsetm
if (i1<0)     {cout<<"err 8353";getch();exit(0);}
if (i2<0)     {cout<<"err 8536";getch();exit(0);}
if (i1>=NFAC) {cout<<"err 7353";getch();exit(0);}
if (i2>=NTRM) {cout<<"err 7536";getch();exit(0);}
if (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)
{
gotoxy(5,14);cout<<"Particle: "<<ch1;
return;
}
void mat(char *ch1)
{
gotoxy(5,16);cout<<"Matrix: "<<ch1;
return;
}

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