//   dsprel6f  Sep 10  small changes
//    dsprel6e   Sep 9  added more materials
//    dsprel6d    Sep 9   try double precision (works); increase speed
//  dsprel6.cpp Aug 31  determine anisotropy of CdS
//   - fixed up random k selector to make it spherically random
//   and also avoid small k.
//   - keep track of directions
//  dsprel5.cpp Aug 31  - test for nu=0.35 and nu=0.15 confirms correct
//     Clong and Cshear for random k directions within 0.1% or so
//     provided k has magnitude in right range.
//  dsprel4   includes k3pt and k6pt, but not checked yet
//  dsprel3  Aug 31
//   2   Aug 31
// disprel1.cpp   Aug 30, 2002  dispersion relation of 3D crystal
//    I am trying to figure out the 3D dispersion relation of a
//    hexagonal lattice.  The five force constants are supplied.
//    In fact I only care about the long wavelength acoustic sound
//    speed.  I just want to know how speed of sound varies with
//    direction and polarization.  (There are three polarizations)

#include <iostream.h>
#include <conio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#include <graphics.h>

double blv[3][3]; // 3 Bravais lattice vectors
double uca,ucc;
double kx,ky,kz,kmag;
int nx,ny,nz;
double dt=0.0001;
double totime=0;
double omega=15;
double u000x,u000y,u000z;
double vel000x,vel000y,vel000z;


/*
//  nu=0.15  (see notes Sun Aug 18 2002 15:58)
//  Use hex2.cpp to get these values
double C11=90e9; // GPa
double C44=37.06e9; // GPa
double densityrho=1000;
double k6pt=-3.287e9;
double ksp1=5.659e10;
double k3pt=-3.909e9;
double ksp3=4.970e10;
double ksp4=2.145e10;
*/

/*
//  nu=0.35  (see notes Sun Aug 18 2002 15:58)
//  Use hex2.cpp to get these values
double C11=90e9; // GPa
double C44=20.77e9; // GPa
double densityrho=1000;
double k6pt=4.298e9;
double ksp1=3.171e10;
double k3pt=5.112e9;
double ksp3=1.690e10;
double ksp4=1.202e10;
*/


/*
// CdSe
double densityrho=5810;
double C11=7.4600e+10;
double C12=4.6100e+10;
double C13=3.9400e+10;
double C33=8.1700e+10;
double C44=1.3000e+10;
double k6pt=4.0973e+09;
double ksp1=2.2535e+10;
double k3pt=9.2628e+09;
double ksp3=2.3245e+10;
double ksp4=7.5253e+09;
*/

/*
// CdS
double densityrho=4825;
double C11=9.0700e+10;
double C12=5.8100e+10;
double C13=5.1000e+10;
double C33=9.3800e+10;
double C44=1.5040e+10;
double k6pt=5.5811e+09;
double ksp1=2.5671e+10;
double k3pt=1.1342e+10;
double ksp3=1.9677e+10;
double ksp4=8.7062e+09;
*/


// Al2O3
double densityrho=3980;
double C11=4.9600e+11;
double C12=1.6400e+11;
double C13=1.1500e+11;
double C33=4.9800e+11;
double C44=1.4800e+11;
double k6pt=-5.1216e+09;
double ksp1=2.6465e+11;
double k3pt=1.8876e+10;
double ksp3=2.6521e+11;
double ksp4=8.5673e+10;


/*
// Wurtzite BN at 300K
double densityrho=3487; // www.ioffe.rssi.ru
double C11=9.8200e+11;
double C12=1.3400e+11;
double C13=7.4000e+10;
double C33=1.0770e+12;
double C44=3.8800e+11;
double k6pt=-4.8733e+10;
double ksp1=6.6977e+11;
double k3pt=-3.8633e+10;
double ksp3=7.5726e+11;
double ksp4=2.2460e+11;
*/
 
/*
// InN 300K
double densityrho=6810; // ioffe.rssi
double C11=1.9000e+11;
double C12=1.0400e+11;
double C13=1.2100e+11;
double C33=1.8200e+11;
double C44=1.0000e+10;
double k6pt=1.7227e+10;
double ksp1=8.6145e+10;
double k3pt=-1.9780e+10;
double ksp3=1.9272e+10;
double ksp4=5.7887e+09;
*/

/*
// InN 300K (Kim '96 calculation Phys Rev B 53 (1996) )
double densityrho=6810; // ioffe.rssi
double C11=2.7100e+11;
double C12=1.2400e+11;
double C13=9.4000e+10;
double C33=2.0000e+11;
double C44=4.6000e+10;
double k6pt=7.4497e+09;
double ksp1=1.2930e+11;
double k3pt=1.0874e+10;
double ksp3=5.4156e+10;
double ksp4=2.6628e+10;
*/

/*
// AlN 300K
double densityrho=3230; // ioffe.rssi  Goldberg 2001;
double C11=4.1000e+11;
double C12=1.4900e+11;
double C13=9.9000e+10;
double C33=3.8900e+11;
double C44=1.2500e+11;
double k6pt=-4.0352e+09;
double ksp1=2.0268e+11;
double k3pt=3.1041e+10;
double ksp3=1.9230e+11;
double ksp4=7.2359e+10;
*/

/*
// Wurtzite GaN 300K
double densityrho=6150;
double C11=3.9000e+11;
double C12=1.4500e+11;
double C13=1.0600e+11;
double C33=3.9800e+11;
double C44=1.0500e+11;
double k6pt=1.5520e+08;
double ksp1=1.9791e+11;
double k3pt=1.7501e+10;
double ksp3=2.0028e+11;
double ksp4=6.0781e+10;
*/

//------------------------------------------------------------
/*
// Wurtzite BN at 300K
double densityrho=3487; // ioffe.rssi
double k6pt=-4.8733e+10;
double ksp1=6.6977e+11;
double k3pt=-3.8633e+10;
double ksp3=7.5726e+11;
double ksp4=2.2460e+11;
*/

/*
// AlN 300K
double densityrho=3230; // ioffe.rssi  Goldberg 2001
// older density measurement was 3225)
double k6pt=-4.0352e+09;
double ksp1=2.0268e+11;
double k3pt=3.1041e+10;
double ksp3=1.9230e+11;
double ksp4=7.2359e+10;
*/

/*
double k6pt=0;
double ksp1=1000;
double k3pt=0;
double ksp3=200;
double ksp4=300;
*/
double fxm[3][3];
double fxn[3][3];
double calcdet(void);
double calcdet2(double lambda);
double srch(double lam1,double lam2,double dlam);
double cellmass,cellvol;
double Clong,Cshear1,Cshear2;
long int maxClong,maxCshear;
double mxclux,mxcluy,mxcluz;
double mnclux,mncluy,mncluz;
double mxcsux,mxcsuy,mxcsuz;
double mncsux,mncsuy,mncsuz;
double sumClong,sumCshear;
int minClong,minCshear;
long int countruns;

void main(void)
{ // main

uca=1.0;
ucc=0.93;
// Define Bravais lattice vectors
blv[0][0]=uca;     blv[0][1]=0;                 blv[0][2]=0;
blv[1][0]=0.5*uca; blv[1][1]=0.5*sqrt(3.0)*uca; blv[1][2]=0;
blv[2][0]=0;       blv[2][1]=0;                 blv[2][2]=ucc;

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

cellvol=uca*uca*sqrt(3.0)*0.5*ucc; // unit cell volume
cellmass=cellvol*densityrho;
double u000mag;
int dispaxis;
double F000x,F000y,F000z;
int ib0,ib1,ib2;
int count[5];

maxClong=0;
maxCshear=0;
sumClong=0;
sumCshear=0;
minClong=15000;
minCshear=15000;

countruns=0;
randomize();
LBagain:
countruns++;
//gotoxy(1,1);cout<<countruns<<" ";
long int lix,liy,liz;
LBchoose:
lix=random(801)-400;
liy=random(801)-400;
liz=random(801)-400;
if ((lix*lix+liy*liy+liz*liz>160000)||(lix*lix+liy*liy+liz*liz<80000)) goto LBchoose;
//kx=0.0003*lix; // Original (f-l-o-a-t)
//ky=0.0003*liy;
//kz=0.0003*liz;
//kx=0.00003*lix;
//ky=0.00003*liy;
//kz=0.00003*liz;
kx=0.00010*lix; // 6e
ky=0.00010*liy;
kz=0.00010*liz;

kmag=sqrt(kx*kx+ky*ky+kz*kz);

u000mag=0.01;
for (dispaxis=0;dispaxis<3;dispaxis++)
{ // dispaxis

u000x=0;
u000y=0;
u000z=0;
if (dispaxis==0) u000x=u000mag;
if (dispaxis==1) u000y=u000mag;
if (dispaxis==2) u000z=u000mag;

// Find total force acting on atom 000:
F000x=0;
F000y=0;
F000z=0;
count[0]=0;count[1]=0;count[2]=0;count[3]=0;count[4]=0;
for (ib0=-1;ib0<=1;ib0++)
for (ib1=-1;ib1<=1;ib1++)
for (ib2=-1;ib2<=1;ib2++)
{ // ib0 ib1 ib2
double r012,r01,r012a,r01a,x1,y1,z1,x1a,y1a,z1a,x0a,y0a,z0a,x2,y2,z2,x0,y0,z0;
x0=0; y0=0; z0=0; // average position of atom 0 0 0
//  (x1 y1 z1) = average position of atom with lattice indices (ib0,ib1,ib2)
x1=ib0*blv[0][0]+ib1*blv[1][0]+ib2*blv[2][0];
y1=ib0*blv[0][1]+ib1*blv[1][1]+ib2*blv[2][1];
z1=ib0*blv[0][2]+ib1*blv[1][2]+ib2*blv[2][2];
//  x1a y1a z1a = actual position of atom
float c1=cos(kx*x1+ky*y1+kz*z1);
x1a=x1+u000x*c1;
y1a=y1+u000y*c1;
z1a=z1+u000z*c1;
x0a=x0+u000x*cos(kx*x0+ky*y0+kz*z0);
y0a=y0+u000y*cos(kx*x0+ky*y0+kz*z0);
z0a=z0+u000z*cos(kx*x0+ky*y0+kz*z0);
r012=(x1-x0)*(x1-x0)+(y1-y0)*(y1-y0)+(z1-z0)*(z1-z0);
r01=sqrt(r012);
r012a=(x1a-x0a)*(x1a-x0a)+(y1a-y0a)*(y1a-y0a)+(z1a-z0a)*(z1a-z0a);
r01a=sqrt(r012a);
double nomdist; double kspn;
int ispr;
for (ispr=1;ispr<=3;ispr++)
{ // ispr
if (ispr==1) {nomdist=uca; kspn=ksp1;}
if (ispr==2) {nomdist=ucc; kspn=ksp3;}
if (ispr==3) {nomdist=sqrt(uca*uca+ucc*ucc); kspn=ksp4;}
if (fabs(r01-nomdist)<.01*uca)
{ // first first neighbour - ksp1
count[ispr]++;
//gotoxy(55,14+ispr);cout<<count[ispr]<<" ksp ";
double Fmag=kspn*(r01a-nomdist);
double Fx,Fy,Fz;
Fx=Fmag*(x1a-x0a)/r01a;
Fy=Fmag*(y1a-y0a)/r01a;
Fz=Fmag*(z1a-z0a)/r01a;
F000x+=Fx;
F000y+=Fy;
F000z+=Fz;
} // first first neighbour - kspn
} // ispr
} // ib0 ib1 ib2

int iway;

// Enumerate all 6-atom clusters:
for (ib0=-1;ib0<=1;ib0++)
for (ib1=-1;ib1<=1;ib1++)
for (ib2=-1;ib2<=1;ib2++)
for (iway=-1;iway<=1;iway=iway+2)
{ // ib0 ib1 ib2 - 6 point force - k6pt
// First, locate center of mass of the cluster:
double cmx,cmy,cmz,xmin,xmax,ymin,ymax,zmin,zmax;
double x1,y1,z1,x2,y2,z2;
cmx=0;
cmy=0;
cmz=0;
xmin=100;xmax=-100;
ymin=100;ymax=-100;
zmin=100;zmax=-100;
int ib0a,ib1a,ib2a,ipt;
for (ib2a=ib2;ib2a<=ib2+1;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
float c1=cos(kx*x1+ky*y1+kz*z1);
x2=x1+u000x*c1;
y2=y1+u000y*c1;
z2=z1+u000z*c1;
cmx=cmx+x2;cmy=cmy+y2;cmz=cmz+z2;
if (x1>xmax) xmax=x1;
if (y1>ymax) ymax=y1;
if (z1>zmax) zmax=z1;
if (x1<xmin) xmin=x1;
if (y1<ymin) ymin=y1;
if (z1<zmin) zmin=z1;
} // ib2a
cmx=cmx/6.0;
cmy=cmy/6.0;
cmz=cmz/6.0;
//putpixel(bx0+sfx*cmx,by0-sfx*cmy,10+iway);

double sumr2=0,r2;
for (ib2a=ib2;ib2a<=ib2+1;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
float c1=cos(kx*x1+ky*y1+kz*z1);
x2=x1+u000x*c1;
y2=y1+u000y*c1;
z2=z1+u000z*c1;
r2=(cmx-x2)*(cmx-x2)+(cmy-y2)*(cmy-y2)+(cmz-z2)*(cmz-z2);
sumr2=sumr2+r2;
} // ib2a, ipt

double Fmag;
double F01x,F01y,F01z,r01;

sumr2=sumr2-6.0*((uca*uca/3.0)+(ucc*ucc/4.0));
Fmag=sumr2*k6pt;

for (ib2a=ib2;ib2a<=ib2+1;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
float c1=cos(kx*x1+ky*y1+kz*z1);
x2=x1+u000x*c1;
y2=y1+u000y*c1;
z2=z1+u000z*c1;
r01=sqrt((x2-cmx)*(x2-cmx)+(y2-cmy)*(y2-cmy)+(z2-cmz)*(z2-cmz));
//F01x=Fmag*(cmx-x2)/r01;
//F01y=Fmag*(cmy-y2)/r01;
//F01z=Fmag*(cmz-z2)/r01;
F01x=Fmag*(cmx-x2); // fixed Aug 18
F01y=Fmag*(cmy-y2); // fixed Aug 18
F01z=Fmag*(cmz-z2); // fixed Aug 18

if ((ib0a==0)&&(ib1a==0)&&(ib2a==0))
{ // force on 000
F000x+=F01x;
F000y+=F01y;
F000z+=F01z;
} // force on 000

} // ib2a, ipt

} // ib0 ib1 ib2 - 6 point force - k6pt

int count333=0;
// Enumerate all 3-atom clusters:
for (ib0=-1;ib0<=1;ib0++)
for (ib1=-1;ib1<=1;ib1++)
for (ib2=-1;ib2<=1;ib2++)
for (iway=-1;iway<=1;iway=iway+2)
{ // ib0 ib1 ib2 - 3 point force - k3pt
// First, locate center of mass of the cluster:
// Aug 18: also locate undistorted max and min of cluster
double cmx,cmy,cmz,xmin,xmax,ymin,ymax,zmin,zmax;
double x1,y1,z1,x2,y2,z2;
cmx=0;
cmy=0;
cmz=0;
xmin=100;xmax=-100;
ymin=100;ymax=-100;
zmin=100;zmax=-100;
int ib0a,ib1a,ib2a,ipt;
for (ib2a=ib2;ib2a<=ib2+0;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
float c1=cos(kx*x1+ky*y1+kz*z1);
x2=x1+u000x*c1;
y2=y1+u000y*c1;
z2=z1+u000z*c1;
cmx=cmx+x2;cmy=cmy+y2;cmz=cmz+z2;
if (x1>xmax) xmax=x1;
if (y1>ymax) ymax=y1;
if (z1>zmax) zmax=z1;
if (x1<xmin) xmin=x1;
if (y1<ymin) ymin=y1;
if (z1<zmin) zmin=z1;
} // ib2a
cmx=cmx/3.0;
cmy=cmy/3.0;
cmz=cmz/3.0;
//putpixel(bx0+sfx*cmx,by0-sfx*cmy,10+iway);

double sumr2=0,r2;
for (ib2a=ib2;ib2a<=ib2+0;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
float c1=cos(kx*x1+ky*y1+kz*z1);
x2=x1+u000x*c1;
y2=y1+u000y*c1;
z2=z1+u000z*c1;
r2=(cmx-x2)*(cmx-x2)+(cmy-y2)*(cmy-y2)+(cmz-z2)*(cmz-z2);
sumr2=sumr2+r2;
} // ib2a, ipt

double Fmag;
double F01x,F01y,F01z,r01;

sumr2=sumr2-3.0*(uca*uca/3.0);
Fmag=sumr2*k3pt;

for (ib2a=ib2;ib2a<=ib2+0;ib2a++)
for (ipt=0;ipt<3;ipt++)
{ // ib2a, ipt
ib0a=ib0; if (ipt==1) ib0a+=iway;
ib1a=ib1; if (ipt==2) ib1a+=iway;
x1=ib0a*blv[0][0]+ib1a*blv[1][0]+ib2a*blv[2][0];
y1=ib0a*blv[0][1]+ib1a*blv[1][1]+ib2a*blv[2][1];
z1=ib0a*blv[0][2]+ib1a*blv[1][2]+ib2a*blv[2][2];
float c1=cos(kx*x1+ky*y1+kz*z1);
x2=x1+u000x*c1;
y2=y1+u000y*c1;
z2=z1+u000z*c1;
r01=sqrt((x2-cmx)*(x2-cmx)+(y2-cmy)*(y2-cmy)+(z2-cmz)*(z2-cmz));
//F01x=Fmag*(cmx-x2)/r01;
//F01y=Fmag*(cmy-y2)/r01;
//F01z=Fmag*(cmz-z2)/r01;
F01x=Fmag*(cmx-x2); // fixed Aug 18
F01y=Fmag*(cmy-y2); // fixed Aug 18
F01z=Fmag*(cmz-z2); // fixed Aug 18
// Aug 18: It should be force on 000 and other inside atoms,
// depending which surface we are looking at:

if ((ib0a==0)&&(ib1a==0)&&(ib2a==0))
{ // force on 000
F000x+=F01x;
F000y+=F01y;
F000z+=F01z;
} // force on 000

} // ib2a, ipt

} // ib0 ib1 ib2 - 3 point force - k3pt



fxm[0][dispaxis]=F000x/u000mag;
fxm[1][dispaxis]=F000y/u000mag;
fxm[2][dispaxis]=F000z/u000mag;

//gotoxy(25+16*dispaxis,5 );cout<<F000x/u000mag<<"  ";
//gotoxy(25+16*dispaxis,6 );cout<<F000y/u000mag<<"  ";
//gotoxy(25+16*dispaxis,7 );cout<<F000z/u000mag<<"  ";

} // dispaxis


double r1,r2,r3,r4,trace,detmat;
detmat=calcdet();  // det of fxm[][]
trace=fxm[0][0]+fxm[1][1]+fxm[2][2];
trace=fabs(trace);
//gotoxy(6,13);cout<<"trace = "<<trace<<" ";
//gotoxy(5,14);cout<<"Det = "<<detmat<<" ";

double lambda,det1,adet1,bestdet=1e35,bestl0,bestl1,bestl2,bestl3,bestl0a;

// Search for lowest eigenvalue, which could not
// exceed 1/3 of trace...
bestl1=srch(0,trace*0.33333,0.025*trace);
bestl1=srch(bestl1-0.05*trace,bestl1+0.05*trace,0.0025*trace);
bestl1=srch(bestl1-0.005*trace,bestl1+0.005*trace,0.00025*trace);
bestl1=srch(bestl1-0.0005*trace,bestl1+0.0005*trace,0.000025*trace);
bestl1=srch(bestl1-0.00005*trace,bestl1+0.00005*trace,0.0000025*trace); // 6e
//gotoxy(6,16);cout<<"bestl(0th): "<<bestl1<<" ";
// Now search for a possible zero below that:
bestl1=srch(0,bestl1-0.05*trace,0.025*trace);
bestl1=srch(bestl1-0.025*trace,bestl1+0.045*trace,0.005*trace);
bestl1=srch(bestl1-0.01*trace,bestl1+0.02*trace,0.002*trace);
bestl1=srch(bestl1-0.02*trace,bestl1+0.02*trace,0.001*trace);
bestl1=srch(bestl1-0.002*trace,bestl1+0.002*trace,0.0001*trace);
bestl1=srch(bestl1-0.0002*trace,bestl1+0.0002*trace,0.00001*trace);
bestl1=srch(bestl1-0.00002*trace,bestl1+0.00002*trace,0.000001*trace); // 6e

// Now search for highest eigenvalue, which
// is bounded below based on known lowest eigenvalue and trace:
// and is also bounded above in a related way:
bestl3=srch((trace-bestl1)*0.5,trace-bestl1,0.025*trace);
bestl3=srch(bestl3-0.05*trace,bestl3+0.05*trace,0.0025*trace);
bestl3=srch(bestl3-0.005*trace,bestl3+0.005*trace,0.00025*trace);
bestl3=srch(bestl3-0.0005*trace,bestl3+0.0005*trace,0.000025*trace); // extra 6e
// Look even higher, just to be sure:
bestl3=srch(bestl3+0.075*trace,trace-bestl1,0.025*trace);
bestl3=srch(bestl3-0.2*trace,bestl3+0.2*trace,0.01*trace);
bestl3=srch(bestl3-0.02*trace,bestl3+0.02*trace,0.001*trace);
bestl3=srch(bestl3-0.002*trace,bestl3+0.002*trace,0.0001*trace);
bestl3=srch(bestl3-0.0002*trace,bestl3+0.0002*trace,0.00001*trace);
bestl3=srch(bestl3-0.00002*trace,bestl3+0.00002*trace,0.000001*trace); // extra 6e

// Now find 2nd eigenvalue:
bestl2=trace-(bestl1+bestl3);
if (bestl2<bestl1)
{ // .. algorithm does not correctly find lowest root 100% of time
double r1=bestl2;
bestl2=bestl1;
bestl1=r1;
} //..
if (bestl3<bestl2)
{ // .. this is probably not needed
double r1=bestl3;
bestl3=bestl2;
bestl2=r1;
} //..

//gotoxy(6,17);cout<<"bestl1: "<<bestl1<<" lowest  ";
//gotoxy(6,18);cout<<"bestl2: "<<bestl2<<" middle  ";
//gotoxy(6,19);cout<<"bestl3: "<<bestl3<<" highest  ";

Cshear1=sqrt(bestl1/cellmass)/kmag;
Cshear2=sqrt(bestl2/cellmass)/kmag;
Clong=sqrt(bestl3/cellmass)/kmag;

sumClong=sumClong+Clong;
sumCshear=sumCshear+0.5*(Cshear1+Cshear2);
if (Clong>maxClong) {maxClong=Clong;mxclux=kx/kmag;mxcluy=ky/kmag;mxcluz=kz/kmag;}
if (Clong<minClong) {minClong=Clong;mnclux=kx/kmag;mncluy=ky/kmag;mncluz=kz/kmag;}
if (Cshear1>maxCshear) {maxCshear=Cshear1;mxcsux=kx/kmag;mxcsuy=ky/kmag;mxcsuz=kz/kmag;}
if (Cshear1<minCshear) {minCshear=Cshear1;mncsux=kx/kmag;mncsuy=ky/kmag;mncsuz=kz/kmag;}
if (Cshear2>maxCshear) {maxCshear=Cshear2;mxcsux=kx/kmag;mxcsuy=ky/kmag;mxcsuz=kz/kmag;}
if (Cshear2<minCshear) {minCshear=Cshear2;mncsux=kx/kmag;mncsuy=ky/kmag;mncsuz=kz/kmag;}

//gotoxy(40,17);cout<<Cshear1<<" m/s   ";
//gotoxy(40,18);cout<<Cshear1<<" m/s   ";
//gotoxy(40,19);cout<<Clong<<" m/s   ";
int i1;

if ((countruns%100)==0)
{ // %5==0
long int li1;
cout.precision(3);
gotoxy(1,1);cout<<countruns<<" ";
gotoxy(20,22);cout<<minCshear<<"    ";
li1=sumCshear/countruns;
gotoxy(35,22);cout<<li1<<"    ";
gotoxy(50,22);cout<<maxCshear<<"    ";
gotoxy(20,24);cout<<minClong<<"    ";
cout.precision(2);
gotoxy(20,25);cout<<mnclux<<" "<<mncluy<<" "<<mncluz<<"   ";
gotoxy(20,23);cout<<mncsux<<" "<<mncsuy<<" "<<mncsuz<<"   ";
cout.precision(3);
li1=sumClong/countruns;
gotoxy(35,24);cout<<li1<<"    ";
gotoxy(50,24);cout<<maxClong<<"    ";
cout.precision(2);
gotoxy(50,25);cout<<mxclux<<" "<<mxcluy<<" "<<mxcluz<<"   ";
gotoxy(50,23);cout<<mxcsux<<" "<<mxcsuy<<" "<<mxcsuz<<"   ";

//gotoxy(40,20);cout<<100*(1-(sqrt(bestl1/cellmass)/kmag)/sqrt(C44/densityrho))<<" %   ";
//gotoxy(40,21);cout<<100*(1-(sqrt(bestl2/cellmass)/kmag)/sqrt(C44/densityrho))<<" %   ";
//gotoxy(40,22);cout<<100*(1-(sqrt(bestl3/cellmass)/kmag)/sqrt(C11/densityrho))<<" %  ";
/*
bestl2=srch(bestl2-.02*trace,bestl2+.02*trace,0.001*trace);
bestl2=srch(bestl2-.002*trace,bestl2+.002*trace,0.0001*trace);
bestl2=srch(bestl2-.0002*trace,bestl2+.0002*trace,0.00001*trace);
gotoxy(46,18);cout<<"bestl2: "<<bestl2<<" middle";
*/

gotoxy(1,21);cout<<"trace/trace = "<<(bestl1+bestl2+bestl3)/trace<<"      ";
//gotoxy(6,22);cout<<"det=? "<<(bestl1*bestl2*bestl3)<<" ";
gotoxy(1,22);cout<<"det/det= "<<(bestl1*bestl2*bestl3)/detmat<<"       ";

gotoxy(1,23);cout<<"Press 'q'";
gotoxy(1,24);cout<<"to exit...";
gotoxy(1,25);cout<<"spacebar to cont.";
} // %5==0

char ch1;
if (!kbhit()) goto LBagain;
ch1=getch();
if (ch1==' ') goto LBagain;

} // main

double calcdet(void)
{ //  calcdet
double ans=0;
ans+=fxm[0][0]*fxm[1][1]*fxm[2][2];
ans+=fxm[1][0]*fxm[2][1]*fxm[0][2];
ans+=fxm[2][0]*fxm[0][1]*fxm[1][2];
ans-=fxm[0][0]*fxm[1][2]*fxm[2][1];
ans-=fxm[1][0]*fxm[2][2]*fxm[0][1];
ans-=fxm[2][0]*fxm[0][2]*fxm[1][1];
return ans;
} //  calcdet

double calcdet2(double lambda)
{ //  calcdet2
int i1,i2;
for (i1=0;i1<3;i1++)
for (i2=0;i2<3;i2++)
{ //.
fxn[i1][i2]=fxm[i1][i2];
} //.
for (i1=0;i1<3;i1++)
fxn[i1][i1]=fxn[i1][i1]+lambda;

double ans=0;
ans+=fxn[0][0]*fxn[1][1]*fxn[2][2];
ans+=fxn[1][0]*fxn[2][1]*fxn[0][2];
ans+=fxn[2][0]*fxn[0][1]*fxn[1][2];
ans-=fxn[0][0]*fxn[1][2]*fxn[2][1];
ans-=fxn[1][0]*fxn[2][2]*fxn[0][1];
ans-=fxn[2][0]*fxn[0][2]*fxn[1][1];
return ans;
} //  calcdet2

double srch(double lam1,double lam2,double dlam)
{ //  srch
double lambda,det1,bestdet,lambdaans,adet1;
bestdet=1e35;
for (lambda=lam1;lambda<lam2;lambda+=dlam)
{ //..
det1=calcdet2(lambda);
adet1=fabs(det1);
if (adet1<bestdet) {bestdet=adet1;lambdaans=lambda;}
} //..
return lambdaans;
} //  srch