#define VERN "scp85b.cpp"
double nu_max=30;
double partrad_nm=0.5*9.8; // part radius (in nanometers)
double matrad_nm=64.0*partrad_nm; // matrix radius (in nanometers)
double nu_start=0.1;
double dnu_min=0.00001; // normally 0.00005
// Except for special points like 0 and pi/2, the
// value of theta5 does not affect A, B, C, D, ...
double theta5=0.25*3.1415926;
int modetype=2; // 1=spheroidal  2=torsional
int orderL=1;
int orderm=0;
int hardwallflag=1;
// scp85b   Apr 3  try hard boundary condition too - hardwallflag
// scp85a   Apr 1  try toroidal modes, orderm=0
//   test (SPH,l=0) for AvgSi-AvgSi -- works perfectly
//   test (TOR,l=2) for AvgSi-AvgSi -- works perfectly
// scp81c.cpp  Mar 29  no changes
//  81p9   testing, no important changes
// scp81p8  Mar 25   fix error in multiplybypow()
//   - improve zero search algorithm
// scp81p7.cpp  Mar 25  - can multiply by kLp etc.
//  - change definitions of A B C D E F following standard literature:
//     A old = A / kLm    SPH LONG 1ST KIND matrix
//     B old = B / kTm    SPH TRAN 1ST KIND matrix
//     C old = C / kLm    SPH LONG 2nd KIND matrix
//     D old = D / kTm    SPH TRAN 2nd KIND matrix
//     E old = E / kLp    SPH LONG 1ST KIND particle
//     F old = F / kTp    SPH TRAN 1ST KIND particle
//   - change normdet definition
//    81p6  - no longer keep 'E' positive.  Makes '.3' file
// scp81p5.cpp   - fix spherical Bessel function of 2nd kind (sign error)
//   - exchanged "B" and "C", following Saviot's choice.
//   - partrad and matrad at top.
//   - show modetype and order L after material parameters
//   - file type 3:  modetype and orderL at left
//   March 23  changed parameter in csinginv and print SEV to output file
//  81p4  eval includes theta, write to file.
// scp81p2  Show that stress is continuous at boundary.
// scp81p  Mar 22  Show rr and r-theta stress near matrix surface.
//     (both seem to go to zero)
//  81n3   Plot r and theta components for eigenfunctions to check
//         for continuity  (seems fine) - detar needs to be small enough
// scp81m  try transposing my[][]
// k  check for 0^0 in pow()
// scp81j  Mar 21  - add csinginv()
//  - I don't understand why E and F are so small.  They should
//    be comparable in magnitude to A and B.
// scp81i.cpp  Mar 19  made 81d.txt and 81e.txt reference files
// scp81h   used to make test reference file 81a.txt
// scp81g  Mar 19  (SPH,l=2)
//    - works okay for: Si around Si; thin outer shell;
//      thick outer shell; small inner heavy particle
//  scp81f.cpp  Seems to be working for case of avg Si around avg Si.
//   eta = 2.64 gives 15.0 cm^-1 and 27 cm^-1.
//  d   why is [4][4] 0 ?
// nm[0][]   sin(theta)
// nm[1][]   cos(theta)
// nm[2][]   sin(kLp r)
// nm[3][]   A
// nm[4][]   cos(kLp r)
// nm[5][]   kLp (longitudinal speed of sound in particle)
// nm[6][]   B
// nm[7][]   r
// nm[8][]   sin(kTp r)
// nm[9][]   C
// nm[10][]  cos(kTp r)
// nm[11][]  kTp (transverse speed of sound in particle)
// nm[12][]  D
//    13  unused
// nm[14][]  sin(kLm r)
// nm[15][]  E
// nm[16][]  cos(kLm r)
// nm[17][]  kLm (longitudinal speed of sound in matrix)
// nm[18][]  F
//    19  unused
// nm[20][]  sin(kTm r)
// nm[21][]  G
// nm[22][]  cos(kTm r)
// nm[23][]  kTm (transverse speed of sound in matrix)
// nm[24][]  H

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

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

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

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

int far bt[NFUNC]; // base term
int far nt[NFUNC]; // number of terms
int far mnt[NFUNC]; // max. number of terms
void dispterm(int it2,int ic,int ir);
void ddr(int func2,int func1); //  Take Deriv of func1 with respect to r
void ddth(int func2,int func1); // Take Deriv of func1 with respect to theta
void ddph(int func2,int func1); // Take Deriv of func1 with respect to phi
void gradient(int ior2,int iot2,int iop2,int func);
void curl(int ior2,int iot2,int iop2,int ior1,int iot1,int iop1);
void divergence(int func2,int ior1,int iot1,int iop1);
void cleanup(int func);
void multiply(int func3,int func2,int func1);
void add(int func3,int func2,int func1);
void subtract(int func3,int func2,int func1);
void legendre(int func1,int ordern);
void sphbessel(int func1,int orderL,int imater,int ipolar,int wavetype);
void copyfunc(int func2,int func1);
void stress(int sigrr,int sigrth,int sigrph,int ur,int uth,int uph,int imater);
void multiplyby(int func,double r1);
void multiplyby(int func,complex cr1);
void multiplyby(int func,char *str1);
void divideby(int func,char *str1);
void multiplybypow(int func,char *str1,int ipow);
void settozero(int func);
int ngetm(int i1,int i2);
void nsetm(int i1,int i2,int i3);
double rho[NMATER];
double css[2][NMATER];
double nu1,dnu1;
double frequencyomegar,frequencyomegai;
complex frequencyomegac;
double kwvr[2][NMATER];
double kwvi[2][NMATER];
void evaluate_func(int func2,int func1,double rad_nm,double theta);
int mosttermseveradd=0;
int mosttermseverddr=0;
int mosttermsevercle=0;
char parstring[40];
char matstring[40];
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 cdet(int Nm2y);
#define Nmy 6
complex y3[Nmy][Nmy];
complex my[Nmy][Nmy];
complex amy[Nmy][Nmy+Nmy]; // aux matrix
int colused[Nmy];
complex singvec[Nmy];
double singeigenval;
FILE *outdatfile;
int outdatflag=0;
char outdatfilename[40];
void csinginv(int Nm2y);

void main(void)
{ // main
int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();
randomize();
gotoxy(1,1);cout<<VERN;
gotoxy(50,4);cout<<"Rp = "<<partrad_nm<<" nm ";
gotoxy(50,5);cout<<"Rm = "<<matrad_nm<<" nm";
gotoxy(50,6);cout<<"Rm/Rp = "<<matrad_nm/partrad_nm<<" ";
gotoxy(50,7);cout<<"Modetype = "<<modetype;
gotoxy(50,8);cout<<"L = "<<orderL;
gotoxy(50,9);cout<<"m = "<<orderm;
gotoxy(50,10);cout<<"nu max = "<<nu_max<<" cm^-1";
gotoxy(50,11);cout<<"dnu min = "<<dnu_min<<" cm^-1";
gotoxy(50,12);cout<<"hard wall flag = "<<hardwallflag;

gotoxy(1,2);cout<<"Find nu A B...F for spheroidal modes (make '3:' file - not normalized)";
gotoxy(1,6);cout<<"( .3 file name extension added automatically";
gotoxy(1,7);cout<<"  '/' to skip  'q' to quit";
gotoxy(1,8);cout<<"  'q' to quit )";
gotoxy(1,4);
cout<<"Data file name? ";
cin>>outdatfilename;
if (outdatfilename[0]=='q') {closegraph();return;}
gotoxy(1,6);cout<<"                                            ";
gotoxy(1,7);cout<<"                        ";
gotoxy(1,8);cout<<"                        ";
outdatflag=1;
if (outdatfilename[0]=='/') outdatflag=0;
if (outdatfilename[0]==0) outdatflag=0;
int i1; char ch1,ch2;
i1=strlen(outdatfilename);
ch1=outdatfilename[i1-4];
ch2=outdatfilename[i1-2];
if ((ch1!='.')&&(ch2!='.'))
{ // add ".3" extension
outdatfilename[i1]='.';
outdatfilename[i1+1]='3';
outdatfilename[i1+2]=0;
} // add extension

if (outdatflag) {gotoxy(1,7);cout<<VERN<<"  Data file: "<<outdatfilename<<" ";}
if (outdatflag) {outdatfile=fopen(outdatfilename,"wb");}
if (outdatflag) {fprintf(outdatfile,"%s%s%c%c","3: ",VERN,13,10);}

int bx0=24,by0=480-40;
double sfx4;

int plotwidth=600;
setcolor(6);
line(bx0,479-301,bx0,by0);
line(bx0,by0,bx0+plotwidth,by0);

double c_nms=1.0e-9*2.9979e8; // speed of light in nm/s
double Pi=4.0*atan(1.0);

//rho[0]=2.329; CLP=9.017; CTP=5.372; par("Avg Si"); // Lucien Saviot's average for Si
//rho[0]=2.329; CLP=8.433; CTP=5.843; par("[100] Si");
//rho[0]=2.329; CLP=9.134; CTP=5.843; par("[110]a Si");
//rho[0]=2.329; CLP=9.134; CTP=4.673; par("[110]b Si");
//rho[0]=2.329; CLP=9.356; CTP=5.093; par("[111] Si");
//rho[0]=4.82;  CLP=4.25;  CTP=1.86;  par("CdS");
//rho[0]=5.81;  CLP=3.57;  CTP=1.54;  par("CdSe");
rho[0]=10.50;    CLP=3.747;  CTP=1.740;  par("Avg Ag"); // direction averaged by dspre7e.cpp
//rho[0]=10;    CLP=3.65;  CTP=1.66;  par("Ag");
//rho[0]=8.908; CLP=5.72;  CTP=3.00;  par("Ni");
//rho[0]=5.33;  CLP=5.25;  CTP=3.35;  par("Ge");
//rho[0]=2.2 ;  CLP=5.95;  CTP=3.76;  par("SiO2");
//rho[0]=10;    CLP=3.00;  CTP=1.732;
//rho[0]=0.125; CLP=1;     CTP=0.3;   par("very soft mat'l.");

rho[1]=3.810 ; CLM=4.610; CTM=2.590; mat("BaO-P2O3"); // Voisin et al Physica B 89, 316 (2002)
//rho[1]=2.2 ; CLM=5.95; CTM=3.76; mat("SiO2");
//rho[1]=2.329; CLM=9.017; CTM=5.372; mat("Avg Si"); // Lucien Saviot's average for Si
//rho[1]=4.82; CLM=4.25; CTM=1.86; mat("CdS"); // CdS
//rho[1]=3.6 ; CLM=3.43; CTM=2.14; mat("GeO2"); // GeO2
//rho[1]=0.125;CLM=1;   CTM=0.3;   mat("very soft mat'l.");
//rho[1]=0.25; CLM=0.5; CTM=0.15;  mat("imag. mat'l.");
//rho[1]=0.5;  CLM=1;   CTM=0.3;   mat("imag. mat'l.");
//rho[1]=1;    CLM=2;   CTM=1.2;   mat("imag. mat'l.");
//rho[1]=2;    CLM=4;   CTM=1.2;   mat("imag. mat'l.");
//rho[1]=10;   CLM=30.; CTM=17.32; mat("nu=0.25 mat'l.");
//rho[1]=2.2;  CLM=5.4; CTM=3.332; mat("borosil glass ???");
// There are many kinds of borosilicate glass.

int imater;
for (imater=0;imater<NMATER;imater++)
{ // 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;
} // imater

sfx4=(plotwidth-bx0-0.5*bx0)/nu_max;

int i_nu;

for (i_nu=2;i_nu<=nu_max;i_nu=i_nu+2)
{ //.. i_nu
setcolor(8);
int ix3=bx0+sfx4*(i_nu);
line(ix3,by0-301,ix3,by0);
} //.. i_nu
for (i_nu=10;i_nu<=nu_max;i_nu=i_nu+10)
{ //.. i_nu
setcolor(2);
int ix3=bx0+sfx4*(i_nu);
line(ix3,by0-301,ix3,by0);
} //.. i_nu

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

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

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

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

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

if (modetype==1) {gotoxy(45,8);cout<<"Spheroidal";}
if (modetype==2) {gotoxy(45,8);cout<<"Torsional";}
gotoxy(45,9);cout<<"order L = "<<orderL<<" ";

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

/*
orderL=0;
//sphbessel(0,orderL,MATRIX,LONGIT,Second_kind);
legendre(0,0);
divideby(0,"kLm");

int ic2,irow2,iobj;
ic2=1; irow2=6;
iobj=0;
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;
*/


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


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

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

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

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

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


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

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

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

if (modetype==2)
{ //.. (TOR)
subtract(20,3,6); // phi continuity
subtract(21,9,12); // r-phi force balance
if (hardwallflag!=1) copyfunc(19,12); // r-phi stress in matrix
if (hardwallflag==1) copyfunc(19,6);
nt[9]=-1; nt[12]=-1;
} //.. (TOR)

if ((modetype==1)&&(orderL>=1))
{ // (SPH,L>=1)
subtract(20,1,4); // r continuity
subtract(21,2,5); // theta continuity
subtract(22,7,10); // rr force balance
if (hardwallflag!=1) copyfunc(18,10); // rr stress in matrix
if (hardwallflag==1) copyfunc(18,4); // rr stress in matrix
nt[7]=-1; nt[10]=-1;
subtract(23,8,11); // r-theta force balance
if (hardwallflag!=1) copyfunc(19,11); // r theta stress in matrix
if (hardwallflag==1) copyfunc(19,5);
nt[8]=-1; nt[11]=-1;
// Phi components 0 in any case for spheroidal.
nt[9]=-1;
nt[12]=-1;
} // SPH

if ((modetype==1)&&(orderL==0))
{ // SPH L=0
subtract(20,1,4); // r continuity
subtract(21,7,10); // rr force balance
if (hardwallflag!=1) copyfunc(19,10); // rr stress in matrix
if (hardwallflag==1) copyfunc(19,4);
nt[7]=-1; nt[10]=-1;
// Phi components 0 in any case for spheroidal.
nt[9]=-1;
nt[12]=-1;
} // SPH L=0

int icolor;
DoAgain:
gotoxy(60,8);cout<<"Partic rad "<<partrad_nm<<" nm ";
gotoxy(60,9);cout<<"Matrix rad "<<matrad_nm<<" nm ";
if (outdatflag)
{ //..
fprintf(outdatfile,"%s%s%c%c","Particle: ",parstring,13,10);
fprintf(outdatfile,"%s%8.2f%s%c%c","Rp  ",partrad_nm," nm",13,10);
fprintf(outdatfile,"%s%8.1f%s%c%c","rho ",rho[0]," kg/m^3",13,10);
fprintf(outdatfile,"%s%8.1f%s%c%c","CLP ",CLP," m/s",13,10);
fprintf(outdatfile,"%s%8.1f%s%c%c","CTP ",CTP," m/s",13,10);
fprintf(outdatfile,"%s%s%c%c","Matrix: ",matstring,13,10);
fprintf(outdatfile,"%s%8.2f%s%c%c","Rm  ",matrad_nm," nm",13,10);
fprintf(outdatfile,"%s%8.1f%s%c%c","rho ",rho[1]," kg/m^3",13,10);
fprintf(outdatfile,"%s%8.1f%s%c%c","CLM ",CLM," m/s",13,10);
fprintf(outdatfile,"%s%8.1f%s%c%c","CTM ",CTM," m/s",13,10);
if ((outdatflag)&&(modetype==1)) {fprintf(outdatfile,"%s%d%c%c","Spheroidal  Order l = ",orderL,13,10);}
if ((outdatflag)&&(modetype==2)) {fprintf(outdatfile,"%s%d%c%c","Torsional   Order l = ",orderL,13,10);}
fprintf(outdatfile,"%s%s%s%c%c",
"Amplitudes A...F are normalized by ",VERN," to 9.99 maximum.",13,10);
if (hardwallflag==1)
{ // hard wall
fprintf(outdatfile,"%s%c%c","Hard wall outer boundary condition.",13,10);
} // hard wall
if (hardwallflag==0)
{ //
fprintf(outdatfile,"%s%c%c","Free surface outer boundary condition.",13,10);
} //
fprintf(outdatfile,"%s%c%c",
"==l==m===nu(cm-1)=====A=========B=========C=========D=========E=========F===",13,10);
} //..
icolor=10+random(5);
double maxdetsofar=0.0,normdet,maxfirstdet=0,normdetlast,normdet2last;
double real_det1,real_det1_last;
long int countpoints=0,countfirst=30;
double lowestnu,lowestdet;
normdetlast=0.0;
normdet2last=0.0;
real_det1_last=0.0;

gotoxy(50,10);cout<<"a) First "<<countfirst<<" ";

LB555:

countpoints++;
if (countpoints>countfirst)
{ //...
nu1=nu1+dnu1;
} //...
if (countpoints==countfirst)
{ // starting value of eta
gotoxy(50,11);cout<<"b) nu = "<<nu_start<<" ";
nu1=nu_start;
dnu1=dnu_min;
} // starting value of eta
else if (countpoints<countfirst)
{ //... first 'countfirst' points
nu1=0.95*nu_max*0.001*random(1000);
} //... first 'countfirst' points

// Angular Frequencies (in rad/ns ?)

frequencyomegar=2.0*Pi*c_nms*100.0*nu1;
frequencyomegai=0;

// Wavevectors are in nm^-1 ?
for (imater=0;imater<NMATER;imater++)
{ //.
kwvr[0][imater]=frequencyomegar/css[0][imater];
kwvr[1][imater]=frequencyomegar/css[1][imater];
kwvi[0][imater]=frequencyomegai/css[0][imater];
kwvi[1][imater]=frequencyomegai/css[1][imater];
} //.

int i89,i89b;

int Nm2y;
if ((modetype==1)&&(orderL>=1)) Nm2y=6;
if ((modetype==1)&&(orderL==0)) Nm2y=3;
if (modetype==2) Nm2y=3;

for (i89=0;i89<Nm2y;i89++)
{ //.
nt[7+i89]=-1;
} //.
// Evaluate boundary condition functions at boundaries:

if ((modetype==1)&&(orderL>=1))
{ // SPH
evaluate_func(0+7,20,partrad_nm,theta5); // r continuity
evaluate_func(1+7,21,partrad_nm,theta5); // th continuity
evaluate_func(2+7,22,partrad_nm,theta5); // rr stress difference
evaluate_func(3+7,23,partrad_nm,theta5); // r th stress
evaluate_func(4+7,18,matrad_nm,theta5);  // rr stress at matrix outer surface
evaluate_func(5+7,19,matrad_nm,theta5);  // r th stress at matrix outer surface
} // SPH

if ((modetype==1)&&(orderL==0))
{ // SPH L=0
evaluate_func(0+7,20,partrad_nm,theta5); // r continuity
evaluate_func(1+7,21,partrad_nm,theta5); // rr stress difference
evaluate_func(2+7,19,matrad_nm,theta5);  // rr stress at matrix outer surface
} // SPH L=0

if (modetype==2)
{ // SPH L=0
evaluate_func(0+7,20,partrad_nm,theta5); //
evaluate_func(1+7,21,partrad_nm,theta5); //
evaluate_func(2+7,19,matrad_nm,theta5);  //
} // SPH L=0

for (i89=0;i89<Nm2y;i89++)
{ //..
cleanup(i89+7);
} //..

for (i89=0;i89<Nm2y;i89++)
{ // i89
my[i89][0]=coefof(i89+7,'A');
my[i89][1]=coefof(i89+7,'B');
my[i89][2]=coefof(i89+7,'C');
if (Nm2y>=4) my[i89][3]=coefof(i89+7,'D');
if (Nm2y>=5) my[i89][4]=coefof(i89+7,'E');
if (Nm2y>=6) my[i89][5]=coefof(i89+7,'F');
} // i89

complex det1;
det1=cdet(Nm2y);
real_det1=real(det1);
normdet=abs(det1)*nu1;
if (normdet>maxdetsofar) maxdetsofar=normdet;

if (countpoints>countfirst)
{ //... countpoints>countfirst
if (normdet<normdetlast)
{ // falling
if (fabs(normdet2last-normdet)>1.5*fabs(normdetlast-normdet))
{
if (fabs(normdetlast-normdet)<0.1*normdet) dnu1=dnu1*1.5;
}
if (fabs(normdetlast-normdet)>0.2*normdet) dnu1=dnu1/2.0;
if (fabs(normdetlast-normdet)>0.3*normdet) dnu1=dnu1/2.0;
} // falling
else if ((normdetlast>normdet2last)&&(normdet>normdet2last))
{ // consistently rising
if (fabs(normdet2last-normdet)>1.5*fabs(normdetlast-normdet))
{ //
if (fabs(normdetlast-normdet)<0.1*normdet) dnu1=dnu1*1.5;
} //
if (fabs(normdetlast-normdet)>0.2*normdet) dnu1=dnu1/1.5;
} // consistently rising
if ((fabs(normdet2last-normdet)<1.6*fabs(normdetlast-normdet))
  ||(fabs(normdet2last-normdet)>3.1*fabs(normdetlast-normdet))
  ||((normdet2last<normdetlast)&&(normdetlast>normdet)))
{ // peaking off - slow down
//putpixel(bx0+0+sfx4*nu1,0+by0-100,12);
dnu1=dnu1/(1.5*1.5);
} // peaking off - slow down
if (dnu1<dnu_min) dnu1=dnu_min;

//gotoxy(50,15);cout<<fabs(normdetlast-normdet)<<"   ";
//gotoxy(50,16);cout<<normdet<<"   ";
//cout.precision(5);
//gotoxy(50,14);cout<<fabs(normdet2last-normdet)/fabs(normdetlast-normdet)<<"     ";
//gotoxy(50,16);cout<<fabs(normdet2last-normdet)/normdet<<"     ";
//gotoxy(50,17);cout<<fabs(normdetlast-normdet)/normdet<<"     ";
//gotoxy(50,18);cout<<"dnu1= "<<dnu1<<"   ";
//gotoxy(50,19);cout<<"nu1= "<<nu1<<"   ";
} //... countpoints>countfirst

if (countpoints<countfirst)
{ //.. (for autoscale display)
if (normdet>maxfirstdet) maxfirstdet=normdet;
} //.. (for autoscale display)

if (countpoints>countfirst)
{ // ..countpoints>countfirst
//putpixel(bx0+0+sfx4*nu1,0+by0-440.0*normdetlast/maxfirstdet,4);
//putpixel(bx0+0+sfx4*nu1,0+by0-440.0*normdet2last/maxfirstdet,1);
putpixel(bx0+0+sfx4*nu1,0+by0-440.0*normdet/maxfirstdet,11);
//if (normdet<normdetlast)
//{ //..
//putpixel(bx0+0+sfx4*nu1,0+by0-440.0*normdet/maxfirstdet,13);
//} //..

//putpixel(bx0+1+sfx4*nu1,0+by0-440.0*normdet/maxfirstdet,8);
//putpixel(bx0+1+sfx4*nu1,1+by0-440.0*normdet/maxfirstdet,8);
//putpixel(bx0+0+sfx4*nu1,1+by0-440.0*normdet/maxfirstdet,8);
if (real_det1_last*real_det1<=0)
{ // det sign changed
putpixel(bx0+sfx4*nu1,by0+10,10);
putpixel(bx0+sfx4*nu1,by0+11,10);
putpixel(bx0+sfx4*nu1,by0+12,10);
putpixel(bx0+sfx4*nu1,by0+13,10);
singvec[3]=0;
singvec[4]=0;
singvec[5]=0;
csinginv(Nm2y);  // Singular matrix inverse
{ //..
cout.precision(2);
int ir,ic;
} //..
cout.precision(4);
gotoxy(1,1);cout<<"nu = "<<nu1<<" cm^-1 ";
gotoxy(1,2 );cout<<"sev = "<<singeigenval<<"      ";
gotoxy(1,3 );cout<<"A = "<<real(singvec[0])<<"    ";
gotoxy(1,4 );cout<<"B = "<<real(singvec[1])<<"    ";
gotoxy(1,5 );cout<<"C = "<<real(singvec[2])<<"    ";
if (Nm2y==6)
{ //
gotoxy(1,6 );cout<<"D = "<<real(singvec[3])<<"    ";
gotoxy(1,7 );cout<<"E = "<<real(singvec[4])<<"    ";
gotoxy(1,8 );cout<<"F = "<<real(singvec[5])<<"    ";
} //
if (outdatflag)
{ //.. write to 3: file
fprintf(outdatfile,
"%1d%2d%3d%10.4f%10.6f%10.6f%10.6f%10.6f%10.6f%10.6f%s%9.2e%c%c",
modetype,orderL,orderm,nu1,
real(singvec[0]),
real(singvec[1]),
real(singvec[2]),
real(singvec[3]),
real(singvec[4]),
real(singvec[5]),
"  SEV= ",singeigenval,13,10);
} //.. write to 3: file
/*
setcolor(7);
line(bx0+sfx4*nu1-10,by0-150,bx0+sfx4*nu1+10,by0-150);
putpixel(bx0+sfx4*nu1,by0-150-100*real(singvec[0]),15);
setcolor(15);circle(bx0+sfx4*nu1,by0-150-100*real(singvec[0]),6);
putpixel(bx0+sfx4*nu1,by0-150-100*real(singvec[1]),14);
setcolor(14);circle(bx0+sfx4*nu1,by0-150-100*real(singvec[1]),5);
putpixel(bx0+sfx4*nu1,by0-150-100*real(singvec[2]),13);
setcolor(13);circle(bx0+sfx4*nu1,by0-150-100*real(singvec[2]),4);
putpixel(bx0+sfx4*nu1,by0-150-100*real(singvec[3]),12);
setcolor(12);circle(bx0+sfx4*nu1,by0-150-100*real(singvec[3]),3);
putpixel(bx0+sfx4*nu1,by0-150-100*real(singvec[4]),11);
setcolor(11);circle(bx0+sfx4*nu1,by0-150-100*real(singvec[4]),2);
putpixel(bx0+sfx4*nu1,by0-150-100*real(singvec[5]),10);
setcolor(10);circle(bx0+sfx4*nu1,by0-150-100*real(singvec[5]),1);
*/
double r_nm,ur,utheta,sfx5=630/matrad_nm,sfy5=50,dr_nm=0.1;
setcolor(8);
line(sfx5*partrad_nm,240-150,sfx5*partrad_nm,240+150);
line(sfx5*matrad_nm,240-150,sfx5*matrad_nm,240+150);
/*
for (r_nm=0.1;r_nm<partrad_nm;r_nm=r_nm+dr_nm)
{ // r_nm

evaluate_func(13,1,r_nm,theta5);
cleanup(13);
complex ce,cf;
ce=coefof(13,'E');
cf=coefof(13,'F');
ur=real(ce)*real(singvec[4])+real(cf)*real(singvec[5]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5*ur,14);
nt[13]=-1;
} // r_nm
*/
/*
for (r_nm=partrad_nm;r_nm<matrad_nm;r_nm=r_nm+dr_nm)
{ // r_nm

evaluate_func(13,4,r_nm,theta5);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
ur=real(ca)*real(singvec[0])
  +real(cb)*real(singvec[1])
  +real(cc)*real(singvec[2])
  +real(cd)*real(singvec[3]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5*ur,14);
nt[13]=-1;
} // r_nm
*/
/*
for (r_nm=0.1;r_nm<partrad_nm;r_nm=r_nm+dr_nm)
{ // r_nm

evaluate_func(13,2,r_nm,theta5);
cleanup(13);
complex ce,cf;
ce=coefof(13,'E');
cf=coefof(13,'F');
utheta=real(ce)*real(singvec[4])+real(cf)*real(singvec[5]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5*utheta,12);
nt[13]=-1;
} // r_nm
*/
/*
for (r_nm=partrad_nm;r_nm<matrad_nm;r_nm=r_nm+dr_nm)
{ // r_nm

evaluate_func(13,5,r_nm,theta5);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
utheta=real(ca)*real(singvec[0])
      +real(cb)*real(singvec[1])
      +real(cc)*real(singvec[2])
      +real(cd)*real(singvec[3]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5*utheta,12);
nt[13]=-1;
} // r_nm
*/
double sigmarr,sigmarth,sfy5s=4e-7;
// Show matrix rr stress:
/*
for (r_nm=0.8*matrad_nm;r_nm<matrad_nm;r_nm=r_nm+0.1)
{ // r_nm

evaluate_func(13,18,r_nm,theta5);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
sigmarr=real(ca)*real(singvec[0])
       +real(cb)*real(singvec[1])
       +real(cc)*real(singvec[2])
       +real(cd)*real(singvec[3]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5s*sigmarr,13);
nt[13]=-1;
} // r_nm
// Show matrix r-theta stress:
for (r_nm=0.8*matrad_nm;r_nm<matrad_nm;r_nm=r_nm+0.1)
{ // r_nm

evaluate_func(13,19,r_nm,theta5);
cleanup(13);
complex ca,cb,cc,cd;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
sigmarth=real(ca)*real(singvec[0])
        +real(cb)*real(singvec[1])
        +real(cc)*real(singvec[2])
        +real(cd)*real(singvec[3]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5s*sigmarth,11);
nt[13]=-1;
} // r_nm
// Show interface stress function:
for (r_nm=0.7*partrad_nm;r_nm<1.3*partrad_nm;r_nm=r_nm+dr_nm)
{ // r_nm

evaluate_func(13,22,r_nm,theta5);
cleanup(13);
complex ca,cb,cc,cd,ce,cf;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
ce=coefof(13,'E');
cf=coefof(13,'F');
sigmarr=real(ca)*real(singvec[0])
       +real(cb)*real(singvec[1])
       +real(cc)*real(singvec[2])
       +real(cd)*real(singvec[3])
       +real(ce)*real(singvec[4])
       +real(cf)*real(singvec[5]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5s*sigmarr,9);
nt[13]=-1;
} // r_nm
for (r_nm=0.7*partrad_nm;r_nm<1.3*partrad_nm;r_nm=r_nm+dr_nm)
{ // r_nm

evaluate_func(13,23,r_nm,theta5);
cleanup(13);
complex ca,cb,cc,cd,ce,cf;
ca=coefof(13,'A');
cb=coefof(13,'B');
cc=coefof(13,'C');
cd=coefof(13,'D');
ce=coefof(13,'E');
cf=coefof(13,'F');
sigmarth=real(ca)*real(singvec[0])
        +real(cb)*real(singvec[1])
        +real(cc)*real(singvec[2])
        +real(cd)*real(singvec[3])
        +real(ce)*real(singvec[4])
        +real(cf)*real(singvec[5]);
putpixel(sfx5*r_nm,240,8);
putpixel(sfx5*r_nm,240+sfy5s*sigmarth,7);
nt[13]=-1;
} // r_nm
*/
//char ch1;
//gotoxy(1,25);cout<<"press key to continue ";ch1=getch();
//if (ch1=='q') return;
//gotoxy(1,25);cout<<"                      ";
} // det sign changed
} // ..countpoints>countfirst


normdet2last=normdetlast;
normdetlast=normdet;
real_det1_last=real_det1;

if ((nu1<nu_max)&&(!kbhit())) goto LB555;

gotoxy(40,25);cout<<"Done...  Press any key to exit.";
if (outdatflag) fclose(outdatfile);
ch1=getch();
return;
} // main

//===========================================================================

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

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

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

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

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

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

return;
} // dispterm

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

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

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

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

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

void cleanup(int func)
{ // cleanup

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

void sphbessel(int func1,int orderL,int imater,int ipolar,int wavetype)
{ // sphbessel
// Mar 24 '03 - fixed sign error for sph bess 2nd kind
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 kind (sign error corrected Mar 24 '03)
cr[it1]=-1; ci[it1]=0; // New line: Mar 24 '03
nsetm(4+6*(kfactor-1),it1,1); // cos(r)
nsetm(5+6*(kfactor-1),it1,-1); // k
} // 1 = 2nd kind
nt[ifaux]=1;
if ((wavetype==2)||(wavetype==3))
{ // 3 = outgoing 2 = ingoing
it1=bt[ifaux]+1;
for (ifac=0;ifac<NFACT;ifac++) nsetm(ifac,it1,0);
cr[it1]=0;
if (wavetype==3) ci[it1]=1; // outgoing
if (wavetype==2) ci[it1]=-1; // ingoing
nsetm(7, it1, -1); // r
nsetm(4+6*(kfactor-1),it1,1); // cos(r)
nsetm(5+6*(kfactor-1),it1,-1); // k
nt[ifaux]=2;
} // 3 = outgoing  2 = ingoing
for (iL=0;iL<orderL;iL++)
{ //. iL
ddr(func1,ifaux);
it1=bt[func1];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //..
nsetm(7,it1,ngetm(7,it1)-1); // r //
nsetm(5+6*(kfactor-1),it1,ngetm(5+6*(kfactor-1),it1)-2); // k factor
} //..
copyfunc(ifaux,func1);
} //. iL
for (iL=0;iL<orderL;iL++)
{ //. iL
it1=bt[ifaux];
for (it1=bt[ifaux];it1<bt[ifaux]+nt[ifaux];it1++)
{ //..
cr[it1]=cr[it1]*(-1); ci[it1]=ci[it1]*(-1);
nsetm(7,it1,ngetm(7,it1)+1); // r //
nsetm(5+6*(kfactor-1),it1,ngetm(5+6*(kfactor-1),it1)+1); // k factor
} //..
} //. iL
copyfunc(func1,ifaux);
nt[ifaux]=-1;
return;
} // sphbessel

void evaluate_func(int func2,int func1,double rad_nm,double theta)
{ // evaluate_func
if (nt[func2]!=-1) {cout<<"* Err evalfunc 0: func "<<func2<<" in use.";getch();exit(0);}
if (mnt[func2]<nt[func1])
 {cout<<"* Err evalfunc 1 mnt<nt";getch();exit(0);}
double r1,r2,r3;
int it1,it2,ifac;
// Step 1: Copy terms from func1 to func2
it1=bt[func1];
nt[func2]=nt[func1];
for (it2=bt[func2];it2<bt[func2]+nt[func2];it2++)
{ //. it2
cr[it2]=cr[it1]; ci[it2]=ci[it1];
int flag3;
for (ifac=0;ifac<NFACT;ifac++)
{ //..
nsetm(ifac,it2,ngetm(ifac,it1));
flag3=1;
// Only factors A, B, C... are kept after evaluation
if ((ifac>=3)&&(((ifac-3)%3)==0)) flag3=0;
if (flag3==1) nsetm(ifac,it2,0);
} //..
it1++;
} //. it2
int ipol,imat,i7;
// Step 2: !!!!!
it2=bt[func2];
for (it1=bt[func1];it1<bt[func1]+nt[func1];it1++)
{ //. it1
if (ngetm(0,it1)!=0)
{ // sin(theta)
cr[it2]=cr[it2]*pow(sin(theta),ngetm(0,it1));
ci[it2]=ci[it2]*pow(sin(theta),ngetm(0,it1));
} // sin(theta)
if (ngetm(1,it1)!=0)
{ // cos(theta)
cr[it2]=cr[it2]*pow(cos(theta),ngetm(1,it1));
ci[it2]=ci[it2]*pow(cos(theta),ngetm(1,it1));
} // cos(theta)
if (ngetm(7,it1)!=0)
{ // r
r1=rad_nm;
if ((r1==0.0)&&(ngetm(7,it1)==0)) {cout<<"* err pow 1234 *";exit(0);}
r2=pow(r1,ngetm(7,it1));
cr[it2]=cr[it2]*r2; ci[it2]=ci[it2]*r2;
} // r
for (ipol=0;ipol<2;ipol++)
{ // ipol
for (imat=0;imat<NMATER;imat++)
{ // imat
i7=6*ipol+12*imat;
if (ngetm(i7+5,it1)!=0)
{ //..
complex cr2,cr3;
cr2=complex(kwvr[ipol][imat],kwvi[ipol][imat]);
if ((cr2==0.0)&&(ngetm(i7+5,it1)==0)) {cout<<"* err pow 1235 *";exit(0);}
cr2=pow(cr2,ngetm(i7+5,it1));
cr3=complex(cr[it2],ci[it2]);
cr3=cr2*cr3;
cr[it2]=real(cr3);
ci[it2]=imag(cr3);
} //..
if (ngetm(i7+2,it1)!=0)
{ //.
complex cr1,cr2,cr3;
cr1=complex(kwvr[ipol][imat]*rad_nm,kwvi[ipol][imat]*rad_nm);
cr2=sin(cr1);
if ((cr2==0.0)&&(ngetm(i7+2,it1)==0)) {cout<<"* err pow 1236 *";exit(0);}
cr2=pow(cr2,ngetm(i7+2,it1));
cr3=complex(cr[it2],ci[it2]);
cr3=cr2*cr3;
cr[it2]=real(cr3);
ci[it2]=imag(cr3);
} //.
if (ngetm(i7+4,it1)!=0)
{ //..
complex cr1,cr2,cr3;
cr1=complex(kwvr[ipol][imat]*rad_nm,kwvi[ipol][imat]*rad_nm);
cr2=cos(cr1);
if ((cr2==0.0)&&(ngetm(i7+4,it1)==0)) {cout<<"* err pow 1237 *";exit(0);}
cr2=pow(cr2,ngetm(i7+4,it1));
cr3=complex(cr[it2],ci[it2]);
cr3=cr2*cr3;
cr[it2]=real(cr3);
ci[it2]=imag(cr3);
} //..
} // imat
} // ipol
it2++;
} //. it1
return;
} // evaluate_func

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

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

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

//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)
{ //
strcpy(parstring,ch1);
//gotoxy(5,15);cout<<"Particle: "<<ch1; return;
} //
void mat(char *ch1)
{ // mat
strcpy(matstring,ch1);
//gotoxy(5,17);cout<<"Matrix: "<<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


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

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

void csinginv(int Nm2y)
{ // csinginv - inverse of matrix, even if singular
char ch1;
int i1,i2,ir,ic,timesthrough,irb,icb;
double maxabs,absr,r1;
complex c1,c2,c3;

for (ir=0;ir<Nm2y;ir++)
for (ic=0;ic<Nm2y;ic++)
{ //. ir ic
amy[ir][ic]=my[ir][ic];
amy[ir][Nm2y+ic]=complex(0.0,0.0);
amy[ir][Nm2y+ir]=complex(1.0,0.0);
} //. ir ic

for (ic=0;ic<Nm2y;ic++) colused[ic]=0;

timesthrough=0;
LB543: // next time through

maxabs=-1.0;
for (ir=0;ir<Nm2y;ir++)
for (ic=0;ic<Nm2y;ic++)
{ //. ir ic
absr=abs(amy[ir][ic]);
if ((absr>maxabs)&&(colused[ic]==0)&&(colused[ir]==0))
{ // ..
irb=ir;
icb=ic;
maxabs=absr;
} // ..
} //. ir ic

for (ic=0;ic<Nm2y+Nm2y;ic++)
{ //. ic
c1=amy[irb][ic];
   amy[irb][ic]=amy[icb][ic];
   amy[icb][ic]=c1;
} //. ic

c2=amy[icb][icb];
//if (abs(c2)<1e-6)
if (abs(c2)<1e-10)
{ //..
//c2=complex(1.0e-6,0.0);
c2=complex(1.0e-10,0.0);
} //..

for (ic=0;ic<Nm2y+Nm2y;ic++)
{ //. ic
amy[icb][ic]=amy[icb][ic]/c2;
} //. ic

for (ir=0;ir<Nm2y;ir++)
{ //.. ir
c1=amy[ir][icb]/amy[icb][icb];
if (ir!=icb)
{ //.
for (ic=0;ic<Nm2y+Nm2y;ic++)
{ //. ic
amy[ir][ic]=amy[ir][ic]-c1*amy[icb][ic];
} //. ic
} //.
} //.. ir

colused[icb]=1;
timesthrough++;
if (timesthrough<Nm2y) goto LB543;

// Determine singular vector:
maxabs=-1.0;
for (ir=0;ir<Nm2y;ir++)
for (ic=0;ic<Nm2y;ic++)
{ //. ir ic
absr=abs(amy[ir][ic+Nm2y]);
if (absr>maxabs) { irb=ir; icb=ic; maxabs=absr; }
} //. ir ic

for (ir=0;ir<Nm2y;ir++)
singvec[ir]=amy[ir][icb+Nm2y];

maxabs=-1.0;
for (ic=0;ic<Nm2y;ic++)
{ // ic
absr=abs(singvec[ic]);
if (absr>maxabs) maxabs=absr;
} // ic

// Normalize singular vector for convenience:
if (maxabs>0)
{ //.
for (ic=0;ic<Nm2y;ic++)
{ // ic
singvec[ic]=9.99*singvec[ic]/maxabs;
} // ic
} //.

// Compute singular eigenvalue:
// (If all has gone well, this should be close to zero)
singeigenval=-99.0;
for (ir=0;ir<Nm2y;ir++)
{ // ir
complex c1;
c1=complex(0,0);
for (ic=0;ic<Nm2y;ic++)
{ // ic
c1=c1+my[ir][ic]*singvec[ic];
} // ic
if (abs(singvec[ir])>0.1)
{ //.
r1=abs(c1)/abs(singvec[ir]);
if (r1>singeigenval) singeigenval=r1;
} //.
} // ir
return;
} // csinginv