int drawanimation=1;
int ipybot=400,ipytop=295;
float Clong,Cshear; // long and shear wave velocity (Jul 24)

float uca=1;  // side length of cubic unit cell (m)

#define VERN "cc3bdmp4.cpp"
#define LennJonesflag 0
#define zview 1
#define makedatfile 0

//   4  Jan 26  make animation 40t.gif for md40.htm
// cc3bdmp3  Jan 26   make 3D pictures for md40.htm
// cc3bdmp2  Jan 25  sphere on hemisphere - figures for md40.htm
// cc3bdmp1  Jan 25  try surface damping, optimize experimentally
//   3  Jan 24   experimenting
// cc3bmod2.cpp  Jan 23, 2003
//    - 2 ExtraRowsForTitle
//    - used "iatomtype" to signify what particles are fixed in position
//    - changed "beta" to "eta"; "alpha" to "xi"
//    - it might be a good idea to plot versus wavenumbers, rather than "eta"
//  Aug. 7, 2002   added Y2O3
// cc3bmod  Aug 5  g/cc, GPa input units, display
//   Aug 4 -- oops! I noticed that I set mvx[] = 0 by mistake
//    (That is why the breathing mode worked while while others didn't.)
//  cc3mod.cpp  Aug 3  8pt
//  cc2mod.cpp  Jul 31  ksp1, ksp2 only (Poisson=.25 if isotropic)
//     (tested and works perfectly for breathing, dog, sph2 modes)
//  Jul 30 2002  cc1mod     suspected error in spring constants
//  Jul 25   fix tick labels, .001 tick marks
//  Jul 24 2002  changed Cshear to 0.619 to 0.6947
//   Jul 22  improve scale
//   Jul 21
//   Jul 20   ang4 ang5 random again (previously both 0)
//   Jul 2   plot infinite R asymptote chart
//   Jul 1 02  - watch 3 atoms, brightness chart versus R
//   June 29   Add legendre poly func and "naive" mode, spherical Bessel
//    Ye toroidal.  Numerical spherical bessel function too.
//   mod2  Jun 28  added legendre poly up to order 4
//  md4mod.cpp   June 26
//  md2mod3.cpp  June 24  show "fr" for cursor
//  2  June 8    excitetype;   watch y motion on 2nd atom
// md2mod.cpp  June 7  eps=0.0174...  improvements
//  612m  June 1   vibrational modes
//  612k.cpp  May 31   tip velocity ratio
// 3d612j2 May 31 fixed energy error
// 3d612j  May 30, but need to invert I; May 31 invert I
//   i  angular momentum; energy ratio; back to float
//  h5   show impulse as function of time
//   h4   try float for position, etc.
// 3d612h.cpp   May 25th
// 3d612g2 looking for problem; contact, liftoff, hop
//   - it is important that dt is small enough
//   - integers for velocity is a problem since there are
//     large forces adding up, or something else wrong.
//   - Also, it runs much slower when I use integers for mv...
// 3d612g   last three neighbours within 255, so use u.s. char
//     vscale mvx now int
//   f   May 24    max optional   up to 11700 atoms
// 3d612e Thurs May 16, May 17   r33 white dot
//  3d612c.cpp  May 15  use dynamic allocation to allow
//   up to 9460 molecules... sphere of radius 13.1
//  3d612b.cpp  May 14  May 15
// 612h.cpp  May 12   model of friction as sinusoidal hard floor
//  612g May 9
// 612f.cpp May 9  found size dependence, apparently that coef of rest
//   approaches 1 in limit of large size, since time scale of
//   collision gets longer and longer when sphere is large.
//  (but of course these are circles and not spheres)
//  612e.cpp   (1st neigh only)
//   - put in dissipative term
//  612d.cpp  May 8
//    - interact with near neighbours only (incl 2nd neigh)
//  Demonstrated coef of rest. near 1 over wide range of speed.
//  (circr=24, dt=0.02)
//  May 7  612c.cpp  - but it occurs to me that the actual
//   equilibrium distance is not the same as the two-body
//   equilibrium distance due to the longer range terms.
//   Therefore, the ball is highly excited when it starts.
//  - dt must be small enough
//  - It is possible to speed up the calculation by constucting
//    a near neighbour table, since the object is solid and
//    keeps the same set of near neighbours unless something
//    extreme happens.
//  - What is the speed of sound?  What is Young's modulus?
//  - What is the density?
//  - Extending this to 3D from 2D seems straightforward.
//  - initial settling and cooling period before starting
//    seems to help a lot.
// 612b.cpp  May 6
// elast2  May 6, 2002   viscoelastic model of dice
//    but this is based on linear theory (small vibrations)
//   so I can't see how to make this work...
// Using Lennard Jones 6-12 potential.
// Equilibrium distance is 1.12246 units.
//   Floor is assumed to be infinitely hard.

#include <iostream.h>
#include <conio.h>
#include <graphics.h>
#include <math.h>
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <dos.h>

int MaxNm,ActNm;
int getmnl(int im1,int in2);
void setmnl(int im1,int in2,int im3);

// Molecule positions:
#define MaxNn 9
unsigned int *m0nl,*m1nl,*m2nl;
unsigned int *m3nl,*m4nl,*m5nl;
unsigned int *m6nl,*m7nl,*m8nl;
int *m9nl;
unsigned char *m13nl,*m11nl;

unsigned char *mcol;
unsigned char *oldmx,*oldmz;
char *oldmy;
unsigned char *iatomtype;

void makescr(void);

// and velocities...
float *mx,*my,*mz;
float *mvx,*mvy,*mvz;
float *max,*may,*maz;

float r122,r121,dx,dy,dz,r124,r126,r127,r1212,r1213,r12,r128,r1214;

float mindist2=-1,maxdist2=-1,mindistever=1000;
float totalenergy,mintotalenergy=2000,maxtotalenergy=-2000;
char beepc=7;

float dissipcoef;
long int vmaxvalever=0;

#define nofr 640
float *fsumr1,*fsumi1;
float *fsumr2,*fsumi2;
float *fsumr3,*fsumi3;
unsigned char *oldfgraf;

float LegenP(int n,float x);
float SphBess1j(int n,float x);
float SeriesSphBess1j(int n,float x);
float rPix(int n,float x,float y,float z);
float rPiy(int n,float x,float y,float z);
float rPiz(int n,float x,float y,float z);
float CrPix(int n,float x,float y,float z);
float CrPiy(int n,float x,float y,float z);
float CrPiz(int n,float x,float y,float z);
float C2rPix(int n,float x,float y,float z);
float C2rPiy(int n,float x,float y,float z);
float C2rPiz(int n,float x,float y,float z);
float ksguess=1,klguess=1;
float THETAx(int n,float x,float y,float z);
float THETAy(int n,float x,float y,float z);
float THETAz(int n,float x,float y,float z);
double factorial(double arg);

void mn(char *ch1);
char matername[50];
float ksp1,ksp2,k8pt;
float densityrho,cellmass;
int n8pt[8];

float dampadjustfac=0.5; // rough experimental optimization Jan 25 03 cc3bdmp1.cpp
int dampadjustcol=1;
float nprad=1,npcmz=1;
void circleb(float rix,float riy,float r1,float r2,int icol);
FILE *scrfout;
long int licountframes=0;

void main(void)
{ // main
clrscr();

char scrfilename[40];
gotoxy(65,8);cout<<"Filename";
//gotoxy(65,9);cout<<"? ";cin>>scrfilename;
strcpy(scrfilename,"40t.scr");
scrfout=fopen(scrfilename,"wb");

float C11,C12,C44,bulkmodulus,shearmodulus,Youngsmodulus,Poissonratio,lambdaLame,r1;
float Y100,Poisson100;
int matertype,isotropic=-1,i1;
LBnewm:
strcpy(matername,"?????");
for (i1=8;i1<=25;i1++)
{ //.
gotoxy(40,i1);cout<<"                                ";
} //.
for (i1=10;i1<=25;i1++)
{ //.
gotoxy(10,i1);cout<<"                                ";
} //.
gotoxy(1,25);cout<<"                                ";
gotoxy(1,1);cout<<"Material Property Input Options:";
gotoxy(1,2);cout<<"--------------------------------";
gotoxy(1,3);cout<<"1. Specify density, Y, Poisson ratio";
gotoxy(1,4);cout<<"2. ";
gotoxy(1,5);cout<<"3. Cubic crystal: (dens., C11, C12, C44)";
gotoxy(1,6);cout<<"4. ";
gotoxy(1,5+2); cout<<"5.  Si   crys 300K";
gotoxy(1,6+2); cout<<"6.  Ge   crys 300K";
gotoxy(1,7+2); cout<<"7.  KCl  crys  4 K";
gotoxy(1,8+2); cout<<"8.  LiCl crys 300K";
gotoxy(1,9+2); cout<<"9.  LiF  crys  0 K";
gotoxy(1,10+2);cout<<"10. GaSb crys 300K";
gotoxy(1,11+2);cout<<"11. MgO  crys 300K";
gotoxy(1,12+2);cout<<"12. NaCl crys 300K";
gotoxy(1,13+2);cout<<"13. GaP  crys 300K";
gotoxy(1,14+2);cout<<"14. GaAs crys 300K";
gotoxy(1,15+2);cout<<"15. CsBr crys 300K";
gotoxy(1,16+2);cout<<"16. Diamond   300K";
gotoxy(1,17+2);cout<<"17. Ag   crys 300K";
gotoxy(1,18+2);cout<<"18. InAs crys 293K";
gotoxy(1,19+2);cout<<"19. Fe   crys 300K";
gotoxy(1,20+2);cout<<"20. Ni   crys 300K";
gotoxy(1,21+2);cout<<"21. Al   crys 300K";
gotoxy(1,22+2);cout<<"22. Cu   crys 300K";
gotoxy(1,23+2);cout<<"23. Au   crys 300K";

gotoxy(21,7);  cout<<"24. RbI  crys 300K";
gotoxy(21,8);  cout<<"25. W    crys  0 K";
gotoxy(21,9);  cout<<"26. Ta   crys  0 K";
gotoxy(21,10); cout<<"27. V    crys  0 K";
gotoxy(21,11); cout<<"28. BaF2 crys 300K";
gotoxy(21,12); cout<<"29. InSb crys 300K";
gotoxy(21,13); cout<<"30. RbBr crys 300K";
gotoxy(21,14); cout<<"31. CsI  crys 300K";
gotoxy(21,15); cout<<"32. Nb   crys 300K";
gotoxy(21,16); cout<<"33. RbCl crys 300K";
gotoxy(21,17); cout<<"34. KF   crys 300K";
gotoxy(21,18); cout<<"35. InP  crys 300K";
gotoxy(21,19); cout<<"36. GaN ZincBlende";
gotoxy(21,20); cout<<"37. Mo   crys     ";
gotoxy(21,21); cout<<"38. Y2O3 crys 293K";
gotoxy(21,22); cout<<"39. - - -         ";
gotoxy(21,23); cout<<"40. - - -         ";
gotoxy(21,24); cout<<"41. - - -         ";
gotoxy(21,25); cout<<"42. - - -         ";


gotoxy(41,7); cout<<"43. - - -  ";//
gotoxy(41,8); cout<<"44. - - -  ";//
gotoxy(41,9); cout<<"45. - - -  ";//
gotoxy(41,10);cout<<"46. Copper, Cu";//
gotoxy(41,11);cout<<"47. Aluminum, Al";
gotoxy(41,12);cout<<"48. Cast Iron ";
gotoxy(41,13);cout<<"49. Lead, Pb";
gotoxy(41,14);cout<<"50. Lucite ";//.
gotoxy(41,15);cout<<"51. Silver, Ag";//..
gotoxy(41,16);cout<<"52. Steel ";//
gotoxy(41,17);cout<<"53. fcc crys ";//
gotoxy(41,18);cout<<"54. Brass ";
gotoxy(41,19);cout<<"55. Europium, Eu";//
gotoxy(41,20);cout<<"56. Alumina Al2O3";
gotoxy(41,21);cout<<"57. 96%SilicaGlass";
gotoxy(41,22);cout<<"58. Fused Quartz";
gotoxy(41,23);cout<<"59. Silcon Carbide";
gotoxy(41,24);cout<<"60. Glass (Pyrex)";
gotoxy(41,25);cout<<"61. - - -        ";

gotoxy(61, 7); cout<<"62. - - -        ";
gotoxy(61, 8); cout<<"63. - - -        ";
gotoxy(61, 9); cout<<"64. - - -        ";
gotoxy(61,10); cout<<"65. - - -        ";
gotoxy(61,11); cout<<"66. - - -        ";
gotoxy(61,12); cout<<"67. - - -        ";
gotoxy(61,13); cout<<"68. Beryllium, Be";
gotoxy(61,14); cout<<"69. Cerium, Ce";
gotoxy(61,15); cout<<"70. Dysprosium,Dy";
gotoxy(61,16); cout<<"71. Erbium, Er";
gotoxy(61,17); cout<<"72. Ruthenium, Ru";
gotoxy(61,18); cout<<"73. Holmium, Ho";
gotoxy(61,19); cout<<"74. Plutonium, Pu";
gotoxy(61,20); cout<<"75. Thulium, Tm";
gotoxy(61,21); cout<<"76. Uranium, U";
gotoxy(61,22); cout<<"77. Yttrium, Y";
gotoxy(61,23); cout<<"78. Ytterbium, Yb";
gotoxy(61,24); cout<<"79. Tungsten, W";
gotoxy(61,25); cout<<"80. Osmium, Os";




gotoxy(50,4);cout<<"Material number?        ";
matertype=5;
//gotoxy(50,4);cout<<"Material number? ";cin>>matertype;
LBnextmat:

for (i1=1;i1<=25;i1++)
{ //.
gotoxy(1,i1);cout<<"                       ";
} //.
for (i1=7;i1<=25;i1++)
{ //..
gotoxy(21,i1);cout<<"                  ";
} //..
for (i1=7;i1<=25;i1++)
{ //..
gotoxy(41,i1);cout<<"                  ";
} //..
for (i1=7;i1<=25;i1++)
{ //..
gotoxy(61,i1);cout<<"                 ";
} //..
gotoxy(1,1);cout<<"                                         ";
gotoxy(1,2);cout<<"                                         ";
gotoxy(1,3);cout<<"                                         ";
gotoxy(1,5);cout<<"                                         ";
if (matertype<1) return;
{ //.. crystals
float rd=-1.00; int im;
im=matertype;
// From Table 22.2 on page 447 of "Solid State Physics" by
// Ashcroft and Mermin  1976
if (im==5)  {C11=  166;C12=64;  C44=80;  rd=2.33; mn("Si crys 300K");} //
//                 166.0   64.0     79.6    2.329 -  ioffe
if (im==6)  {C11=  129;C12=48;  C44=67;  rd=5.32; mn("Ge crys 300K");} //
//                 126     44.0     67.7    -  ioffe
// From Table 2 and Table 3 on pages 149 and 150 of C. Kittel
// Introduction to Solid State Physics  4th Ed. 1971
if (im==7)  {C11= 48.3;C12=5.4; C44=6.6; rd=2.038;mn("KCl crys @ 4K");} //
if (im==8)  {C11= 49.4;C12=22.8;C44=24.6;rd=2.07; mn("LiCl crys 300K");} //
if (im==9)  {C11=124.6;C12=42.4;C44=64.9;rd=2.646;mn("LiF crys @ 0K");}
if (im==10) {C11= 88.5;C12=40.4;C44=43.3;rd=5.614;mn("GaSb crys 300K");} //
//                88.3     40.2     43.2    5.614  300K ioffe
if (im==11) {C11=286  ;C12=87  ;C44=148 ;rd=3.60; mn("MgO crys 300K");} //
if (im==12) {C11= 48.7;C12=12.4;C44=12.6;rd=2.17; mn("NaCl crys 300K");} //
if (im==13) {C11=140.5;C12=62.0;C44=70.3;rd=4.138;mn("GaP crys 300K");} //
if (im==14) {C11=119.0;C12=53.4;C44=59.6;rd=5.317;mn("GaAs crys 300K");} //
if (im==15) {C11=30.0; C12=7.8 ;C44=7.6 ;rd=4.43; mn("CsBr crys 300K");} // Kittel Table 3
if (im==16){C11=1076;  C12=125; C44=576; rd=3.516;mn("Diamond 300K");} // Kittel Table 3
//              1079       124      578     3.515   300K ioffe
if (im==17) {C11=124;  C12= 93; C44=46;  rd=10.6; mn("Ag crys 300K");} //Ashcroft& Mermin Table 22.2 - J. DeLaunay 1956
if (im==18) {C11=83.4; C12=45.4;C44=39.5;rd=5.68; mn("InAs crys 293K");}// www.ioffe.rssi.ru - Burenkov et al. 1975
if (im==19) {C11=234;  C12=136; C44=118; rd=7.87; mn("Fe crys 300K");} //Ashcroft& Mermin Table 22.2 - Rayne 1961
if (im==20) {C11=245;  C12=140; C44=125; rd=8.97; mn("Ni crys 300K");} //Ashcroft& Mermin Table 22.2 - Alers 1960
if (im==21) {C11=107;  C12= 61; C44= 28; rd=2.73; mn("Al crys 300K");} // Ashcroft& Mermin Table 22.2 - Ho 1969
if (im==22) {C11=168;  C12=121; C44=75;  rd=9.02; mn("Cu crys 300K");} // Ashcroft& Mermin Table 22.2 - J. DeLaunay 1956
if (im==23) {C11=186;  C12=157; C44=42;  rd=19.5; mn("Au crys 300K");} // Ashcroft/Mermin Table 22.2 - J. DeLaunay 1956
// (im==23) {C11=192.3;C12=163.1;C44=42.0;rd=19.5; mn("Au crys");} // Kittel Table 2
if (im==24) {C11=25.6; C12=3.1; C44=2.9; rd=3.55; mn("RbI cry. 300K");} // Kittel Table 3
if (im==25) {C11=532.6; C12=204.9; C44=163.1; rd=19.317; mn("W crys O K");}
if (im==26) {C11=266.3; C12=158.2; C44=87.5; rd=16.696; mn("Ta crys 0 K");}
if (im==27) {C11=232.4; C12=119.4; C44=46.0; rd=6.051; mn("V crys 0 K");}
if (im==28) {C11=89.1; C12=40.0; C44=25.4; rd=4.886; mn("BaF2 crys 300K");}
if (im==29) {C11=66.7; C12=36.5; C44=30.2; rd=5.77; mn("InSb crys 300K");} // www.ioffe.rssi.ru
// (im==29) {C11=67.2; C12=36.7; C44=30.2; rd=5.77; mn("InSb crys 300K");} // Kittel
if (im==30) {C11=31.7; C12=4.2; C44=3.9; rd=3.35; mn("RbBr crys 300K");}
if (im==31) {C11=24.6; C12=6.7; C44=6.2; rd=4.51; mn("CsI crys 300K");}
if (im==32) {C11=247;  C12=135; C44=28.7; rd=8.58; mn("Nb crys 300K");}
if (im==33) {C11=36.1; C12=6.2; C44=4.7; rd=2.76; mn("RbCl crys 300K");}
if (im==34) {C11=65.6; C12=14.6; C44=12.5; rd=2.48; mn("KF crys 300K");}
if (im==35) {C11=101.1;C12=56.1; C44=45.6; rd=4.81; mn("InP crys 300K");} // www.ioffe.rssi.ru
if (im==36) {C11=293;  C12=159;  C44=155;  rd=6.15; mn("GaN ZincBlende");} // @ 300K - www.ioffe.rssi.ru
//    GaN can also have wurtzite structure, which is hexagonal, not cubic
if (im==37) {C11=46.30; C12=16.10; C44=10.90; rd=10.22; mn("Mo crys");} //
if (im==38) {C11=223.7; C12=112.4; C44=74.6; rd=5.01; mn("Y2O3 crys 293K");} // ref [22] in matprop.htm - Palko et al 2001
if (im==39) {C11=-1; C12=-1; C44=-1; rd=-1; mn("???");}
if (im==40) {C11=-1; C12=-1; C44=-1; rd=-1; mn("???");}
if (im==41) {C11=-1; C12=-1; C44=-1; rd=-1; mn("???");}
if (im==42) {C11=-1; C12=-1; C44=-1; rd=-1; mn("???");}
//  AlxGa1-xAs  has C12<C44 usually- ioffe
//  GaAs1-xSbx  has C12<C44 usually- ioffe

if ((im>=5)&&(im<=42))
{ //.
if (rd<0) goto LBnewm;
C11=C11*1.0e9;
C12=C12*1.0e9;
C44=C44*1.0e9;
densityrho=rd*1000.0;
} //.
} //.. crystals  #5-42

{ // rho Y G  #43-80
float Y=-1,pr=-1,rd=-1; int im;
im=matertype;
if (im==43){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==44){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==45){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==46){rd=8.93; Y=122 ;pr=0.33 ;mn("Copper, Cu");} //
if (im==47){rd=2.7;  Y=71;pr=0.33;mn("Aluminum, Al");} //
if (im==48){rd=7.70; Y=105 ;pr=0.28;mn("Cast iron");} //
if (im==49){rd=11.34;Y=16.5 ;pr=0.43 ;mn("Lead, Pb");} //
if (im==50){rd=1.20; Y=4 ;pr=0.4 ;mn("Lucite");} //
if (im==51){rd=10.5; Y=78 ;pr=0.37 ;mn("Silver, Ag");} //
if (im==52){rd=7.8;  Y=210 ;pr=0.29 ;mn("Steel");} //
if (im==53){rd=0.001;Y=1.215e-9 ;pr=0.2589 ;mn("avg fcc md crystal");} //
if (im==54){rd=8.4;  Y=104 ;pr=0.37 ;mn("Brass");} //
if (im==55){rd=5.26;Y=18.2;pr=0.152;mn("Europium, Eu");}
if (im==56){rd=3.9 ; Y=340 ;pr=0.22 ;mn("Alumina Al2O3");} //
if (im==57){rd=2.18 ; Y=68 ;pr=0.19 ;mn("96% silica glass");} //
if (im==58){rd=2.2 ; Y=70 ;pr=0.17 ;mn("Fused quartz");} //
if (im==59){rd=3.1 ; Y=410 ;pr=0.24 ;mn("Silicon carbide");} // , SiC
if (im==60){rd=2.30; Y=62 ;pr=0.24 ;mn("Pyrex glass");} //
if (im==61){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==62){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==63){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==64){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==65){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==66){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==67){rd=-1 ; Y=-1 ;pr=-1 ;mn("??");}
if (im==68){rd=1.844; Y=303;  pr=0.1667;mn("Beryllium, Be");} //  (nu 0.07-0.18)
if (im==69){rd=6.7 ; Y=41.4;  pr=0.248; mn("Cerium, Ce");} //
if (im==70){rd=8.54 ; Y=61.4; pr=0.237; mn("Dysprosium, Dy");} //
if (im==71){rd=9.05 ; Y=69.9; pr=0.237; mn("Erbium, Er");} //
if (im==72){rd=12.2; Y=414;   pr=0.25;  mn("Ruthenium, Ru");}  //
if (im==73){rd=8.8 ; Y=64.8;  pr=0.231; mn("Holmium, Ho");} //
if (im==74){rd=19.81; Y=96.5; pr=0.18;  mn("Plutonium, Pu");} //  (nu 0.15-0.21)
if (im==75){rd=9.33 ; Y=74;   pr=0.213; mn("Thulium");} //
if (im==76){rd=19.07; Y=190;  pr=0.22;  mn("Uranium, U");} //
if (im==77){rd=4.472; Y=63.5; pr=0.243; mn("Yttrium, Y");} //
if (im==78){rd=6.98 ; Y=23.7; pr=0.207; mn("Ytterbium, Yb");} //
if (im==79){rd=19.3;  Y=345;  pr=0.17;  mn("Tungsten, W");} //
if (im==80){rd=22.6;  Y=552;  pr=0.25;  mn("Osmium, Os");} //

if ((im>=43)&&(im<=100))
{ // get C11,C12,C44
if (rd<0) goto LBnewm;
densityrho=rd*1000.0;
Y=Y*1e9;
C44=Y/(2.0*(1.0+pr));
lambdaLame=Y*pr/((1+pr)*(1-2*pr));
C11=lambdaLame+2*C44;
C12=C11-2*C44; // Assume isotropic
} // get C11,C12,C44
} // rho Y G

if (matertype==1)
{ // 1   Y and Poisson
gotoxy(50-2,11);cout<<"* Density ? (g/cc)         ";
//densityrho=1e3;
gotoxy(50-2,11);cout<<"* Density ? (g/cc) ";cin>>densityrho;
densityrho=densityrho*1000.0;
if (densityrho<=0) {cout<<"* Non-positive density...*";return;}
gotoxy(50-2,11);cout<<"  ";
gotoxy(50-2,15);cout<<"* Young's modulus ? (GPa)        ";
//Youngsmodulus=1e11;
gotoxy(50-2,15);cout<<"* Young's modulus ? (GPa) ";cin>>Youngsmodulus;
Youngsmodulus=Youngsmodulus*1e9;
if (Youngsmodulus<0) {cout<<"* negative Y *";return;}
gotoxy(50-2,15);cout<<"  ";
LBnewPR:
gotoxy(50-2,16);cout<<"* Poisson ratio ?           ";
//Poissonratio=0.15;
gotoxy(50-2,16);cout<<"* Poisson ratio ? ";cin>>Poissonratio;
gotoxy(50-2,16);cout<<"  ";
C44 = Youngsmodulus/(2.0*(1.0+Poissonratio));
lambdaLame=Youngsmodulus*Poissonratio/((1+Poissonratio)*(1-2*Poissonratio));
C11 = lambdaLame + 2*C44;
C12 = C11 - 2*C44; // Assume isotropic
} // 1   Y and Poisson
if (matertype==3)
{ // 3 Cubic crystal
gotoxy(50,11);cout<<"Density ? (g/cc) ";cin>>densityrho;
if (densityrho<=0) {cout<<"* Non-positive density...*";return;}
densityrho=densityrho*1000.0;
gotoxy(50,12);cout<<"C11 ? (GPa) ";cin>>C11;
gotoxy(50,13);cout<<"C12 = lambda (GPa) ? ";cin>>C12;
gotoxy(50,14);cout<<"C44 = G = æ ? (GPa) ";cin>>C44;
C11=C11*1e9;
C12=C12*1e9;
C44=C44*1e9;
} // 3 Cubic crystal
gotoxy(50,9);cout<<matername;
gotoxy(50,10);cout<<"---------------";
cout.precision(4);
gotoxy(50,11);
if (densityrho>800)
 cout<<"Density = "<<0.001*densityrho<<" g/cc   ";
else
 cout<<"Density = "<<densityrho<<" kg/m^3   ";
gotoxy(50,12);
if (C11>1e6)
 cout<<"C11 = "<<1.0e-9*C11<<" GPa  ";
else
 cout<<"C11 = "<<C11<<" Pa  ";
gotoxy(50,13);
if (C12>1e6)
 cout<<"C12 = "<<1.0e-9*C12<<" GPa  ";
else
 cout<<"C12 = "<<C12<<" Pa  ";
gotoxy(50,14);
if (C44>1e6)
 cout<<"C44 = "<<1.0e-9*C44<<" GPa  ";
else
 cout<<"C44 = "<<C44<<" Pa  ";
//ksp1 = uca*(C11-C44);
//ksp2 = uca*(C44-C12);
//if (ksp2<0) {ksp2=fabs(ksp2);gotoxy(65,18);cout<<"ksp2 made pos";}
//if (ksp2<0) {ksp2=0;gotoxy(65,18);cout<<"ksp2 made zero";}

isotropic=0;
if (fabs(C12-(C11-2*C44))<0.00001*fabs(C12)) isotropic=1;
ksp1=uca*(C11-(C12+C44));
ksp2=uca*C44;
k8pt=1.5*uca*(C12-C44);
gotoxy(50,17);cout<<"ksp1 = "<<ksp1<<" N/m    ";
gotoxy(50,18);cout<<"ksp2 = "<<ksp2<<" N/m    ";
gotoxy(50,19);cout<<"k8pt = "<<k8pt<<" N/m    ";
bulkmodulus = (C11+2.0*C12)/3.0;
//shearmodulus = C44; // What I used to think!
shearmodulus = (C11-C12)/2; // correct definition, general to cubic crys. see ioffe.rssi.ru
// From Kittel page 154, problem 4: "Effective Shear Constant"
// "Show that the shear constant (C11-C12)/2 in a cubic crystal is
//  defined by setting exx = -eyy = e/2 ..."
Youngsmodulus = 9.0*bulkmodulus*shearmodulus/(3.0*bulkmodulus+shearmodulus);
Poissonratio = 0.5 * (3.0*bulkmodulus-2.0*shearmodulus)/(3.0*bulkmodulus+shearmodulus);
Y100 = (C11+2*C12)*(C11-C12)/(C11+C12);
Poisson100 = C12/(C11+C12);
// (for Y100 and Poisson 100 see Landau & Lifshitz Theory of Elasticity
// 3rd Edition page 37, solution to problem 3)
lambdaLame=Youngsmodulus*Poissonratio/((1-2*Poissonratio)*(1+Poissonratio));

if (isotropic==1)
{ // isotropic
gotoxy(50,15);
if (Youngsmodulus>1e6)
{ //.
 cout<<"Young mod. = ";
 printf("%5.1f",1.0e-9*Youngsmodulus);
 cout<<" GPa       ";
} //.
else
 cout<<"Young mod. = "<<Youngsmodulus<<" Pa       ";
cout.precision(3);
gotoxy(50,16);cout<<"Poisson ratio = "<<Poissonratio<<"      ";
} // isotropic
if (isotropic==0)
{ // not isotropic
gotoxy(50,15);cout<<"Y[100] = ";
printf("%5.1f",1.0e-9*Y100);
cout<<" GPa       ";
cout.precision(3);
gotoxy(50,16);cout<<"Poisson[100] = "<<Poisson100<<"      ";
} // not isotropic

gotoxy(50,20);
if (bulkmodulus>1e6)
{ //.
 cout<<"B = ";
 printf("%5.1f",1.0e-9*bulkmodulus);
 cout<<" GPa ";
} //.
else
 cout<<"B = "<<bulkmodulus<<" Pa ";
gotoxy(50,21);
if (shearmodulus>1e6)
 cout<<"G, æ, C' = "<<1.0e-9*shearmodulus<<" GPa ";
else
 cout<<"G, æ, C' = "<<shearmodulus<<" Pa ";
gotoxy(50,22);
if (lambdaLame>1e6)
 cout<<"Lambda = "<<1.0e-9*lambdaLame<<" GPa ";
else
 cout<<"Lambda = "<<lambdaLame<<" Pa ";
gotoxy(1,14);cout<<"                 ";
gotoxy(1,15);cout<<"Sound speeds:   ";
gotoxy(1,16);cout<<"-------------";
cout.precision(3);
if (isotropic==0)
{ // not
cout.precision(3);
gotoxy(1,5);cout<<"Zener Anisotropy Ratio:";
gotoxy(1,6);cout<<"A = 2 C44/(C11-C12) = "<<2*C44/(C11-C12)<<" ";
int iCL100=sqrt(C11/densityrho);
int iCT100=sqrt(C44/densityrho);
gotoxy(1,17);cout<<"[100] CL  "<<iCL100<<" m/s ";
gotoxy(1,18);cout<<"[100] CS  "<<iCT100<<" m/s ";
int iCL110=sqrt((0.5*(C11+C12+2.0*C44))/densityrho);
gotoxy(1,19);cout<<"[110] CL  "<<iCL110<<" m/s ";
int iCT110a=sqrt((C44)/densityrho);
gotoxy(1,20);cout<<"[110] CS1 "<<iCT110a<<" m/s ";
int iCT110b=sqrt((0.5*(C11-C12))/densityrho);
gotoxy(1,21);cout<<"[110] CS2 "<<iCT110b<<" m/s ";
int iCL111=sqrt( ((C11+2*C12+4*C44)/3.0) /densityrho);
gotoxy(1,22);cout<<"[111] CL  "<<iCL111<<" m/s ";
int iCT111=sqrt( ((C11-C12+C44)/3.0) /densityrho);
gotoxy(1,23);cout<<"[111] CS  "<<iCT111<<" m/s ";
} // not isotropic
if (isotropic==1)
{ // isotropic
int iCL100=sqrt(C11/densityrho); float CL100=iCL100;
int iCT100=sqrt(C44/densityrho); float CT100=iCT100;
if (iCL100<100) CL100=sqrt(C11/densityrho);
if (iCT100<100) CT100=sqrt(C44/densityrho);
gotoxy(1,17);cout<<"C longi  "<<CL100<<" m/s ";
gotoxy(1,18);cout<<"C shear  "<<CT100<<" m/s ";
} // isotropic

/*
gotoxy(1,25);cout<<"Press spacebar to continue...";
{
char ch1=getch();
if (ch1=='m') {clrscr();goto LBnewm;}
if (ch1=='n') {matertype=matertype+1;clrscr();goto LBnextmat;}
if (ch1=='p') goto LBnewPR;
if (ch1=='q') {closegraph();return;}
}
*/

int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();
struct palettetype pal;
getpalette(&pal);
int ic;
for (ic=0;ic<16;ic++)
{ //
setrgbpalette(pal.colors[ic],4*ic,4*ic,4*ic);
} //

setrgbpalette(pal.colors[7],40,40,63); // for text
setrgbpalette(pal.colors[1],63,40,40);
setrgbpalette(pal.colors[0],63,63,47); // pale yellow ffffc0
gotoxy(1,1);cout<<VERN;
gotoxy(1,2);cout<<"Normal mode frequency from Fourier analysis ";

fsumr1=new float[nofr];
fsumi1=new float[nofr];
fsumr2=new float[nofr];
fsumi2=new float[nofr];
fsumr3=new float[nofr];
fsumi3=new float[nofr];
oldfgraf=new unsigned char[nofr];

FILE *fout; char foutname[30];

if (makedatfile)
{ //
gotoxy(1,10);cout<<"File name? ";cin>>foutname;
fout=fopen(foutname,"wb");
} //

dissipcoef=0;

gotoxy(55,1);cout<<matername<<" ";

//gotoxy(3,5);cout<<"Space for 13200 inside C++ IDE ";
//gotoxy(3,6);cout<<"Space for 15200 when run from DOS ";
MaxNm=7500;
//gotoxy(3,3);cout<<"Max no. of atoms? ";cin>>MaxNm;
if (MaxNm<=0) {closegraph();return;}

//if (!(mnn=new unsigned char[MaxNm])) {cout<<"Not enough for mnn";getch();closegraph();return;}
if (drawanimation==1)
{ //.
if (!(mcol=new unsigned char[MaxNm])) {cout<<"Not enough for mcol";getch();closegraph();return;}
if (!(oldmx=new unsigned char[MaxNm])) {cout<<"Not enough for oldmx";getch();closegraph();return;}
if (!(oldmy=new unsigned char[MaxNm])) {cout<<"Not enough for oldmy";getch();closegraph();return;}
if (!(oldmz=new unsigned char[MaxNm])) {cout<<"Not enough for oldmz";getch();closegraph();return;}
} //
if (!(iatomtype=new unsigned char[MaxNm])) {cout<<"Not enough for iatomtype";getch();closegraph();return;}
if (!(m0nl=new unsigned int[MaxNm])) {cout<<"Not enough for m0nl";getch();closegraph();return;}
if (!(m1nl=new unsigned int[MaxNm])) {cout<<"Not enough for m1nl";getch();closegraph();return;}
if (!(m2nl=new unsigned int[MaxNm])) {cout<<"Not enough for m2nl";getch();closegraph();return;}
if (!(m3nl=new unsigned int[MaxNm])) {cout<<"Not enough for m3nl";getch();closegraph();return;}
if (!(m4nl=new unsigned int[MaxNm])) {cout<<"Not enough for m4nl";getch();closegraph();return;}
if (!(m5nl=new unsigned int[MaxNm])) {cout<<"Not enough for m5nl";getch();closegraph();return;}
if (!(m6nl=new unsigned int[MaxNm])) {cout<<"Not enough for m5nl";getch();closegraph();return;}
if (!(m7nl=new unsigned int[MaxNm])) {cout<<"Not enough for m5nl";getch();closegraph();return;}
if (!(m8nl=new unsigned int[MaxNm])) {cout<<"Not enough for m5nl";getch();closegraph();return;}
if (!(m9nl=new int[MaxNm])) {cout<<"Not enough for m9nl";getch();closegraph();return;}
if (!(m11nl=new unsigned char[MaxNm])) {cout<<"Not enough for m11nl";getch();closegraph();return;}
if (!(m13nl=new unsigned char[MaxNm])) {cout<<"Not enough for m13nl";getch();closegraph();return;}

if (!(mx=new float[MaxNm])) {cout<<"Not enough for mx";getch();closegraph();return;}
if (!(my=new float[MaxNm])) {cout<<"Not enough for my";getch();closegraph();return;}
if (!(mz=new float[MaxNm])) {cout<<"Not enough for mz";getch();closegraph();return;}
if (!(mvx=new float[MaxNm])) {cout<<"Not enough for mvx";getch();closegraph();return;}
if (!(mvy=new float[MaxNm])) {cout<<"Not enough for mvy";getch();closegraph();return;}
if (!(mvz=new float[MaxNm])) {cout<<"Not enough for mvz";getch();closegraph();return;}
if (dissipcoef!=0)
{ //..
if (!(max=new float[MaxNm])) {cout<<"Not enough for max";getch();closegraph();return;}
if (!(may=new float[MaxNm])) {cout<<"Not enough for may";getch();closegraph();return;}
if (!(maz=new float[MaxNm])) {cout<<"Not enough for maz";getch();closegraph();return;}
} //..
gotoxy(1,1);cout<<"Memory space is sufficient... ";

int bx0=370,by0=100;

float dt,ttot;
dt=0.6*uca/sqrt(C11/densityrho);
cout.precision(3);

Clong=sqrt(C11/densityrho);
Cshear=sqrt(C44/densityrho);

float ang4,ang5;

randomize();

float xmax=68;
float ymax=xmax*480/640.0;
float sfx=(640.0-100)/xmax;
int ix,iy,iz,im,im1,im2;

float ObjRad,maxpossrad;
float aspectr,radratio;

int whichshape;
gotoxy(1,5);cout<<"0=sphere";
gotoxy(1,6);cout<<"1=hollow sphere";
gotoxy(1,7);cout<<"2=cylinder";
gotoxy(1,8);cout<<"3=flat cube";
gotoxy(1,9);cout<<"4=sphere on half sphere";
gotoxy(55,2);cout<<"What shape? ";
whichshape=4;
//cin>>whichshape;
if (whichshape<0) {closegraph();return;}
aspectr=1; radratio=1;
if (whichshape==1) {gotoxy(20,6);cout<<"rad ratio? (rin/rout) ";cin>>radratio;}
if (whichshape==2) {gotoxy(15,7);cout<<"h/D ? ";cin>>aspectr;}
if (whichshape==3) {gotoxy(15,8);cout<<"h/w ? ";cin>>aspectr;}
if (whichshape==4)
{ // sphere stuck to half sphere
nprad=0.5;
//gotoxy(15,9); cout<<"radius:   nprad ? ";cin>>nprad;
npcmz=0.45;
//gotoxy(15,10);cout<<"cm shift: npcmz ? ";cin>>npcmz;
} // sphere stuck to half sphere


float cmvx,cmvy,cmvz,cmx,cmy,cmz;
if (whichshape==0) maxpossrad=pow(MaxNm/(1.333333*3.1415926),0.333333);
if (whichshape==4) maxpossrad=pow((1.8*MaxNm)/(1.333333*3.1415926),0.333333); // half sphere
if (whichshape==1) maxpossrad=pow((MaxNm/(1.0-pow(radratio,3.0)))/(1.333333*3.1415926),0.3333333);
//  MaxNm = 3.14159 * maxpossrad^3 * 2.0* aspectr
if (whichshape==2) maxpossrad=pow(MaxNm/(2.0*3.1415926*aspectr),0.333333);
cout.precision(2);
gotoxy(40,7);cout<<"Max. poss.";
gotoxy(40,8);cout<<"rad is "<<maxpossrad;
ObjRad=5.0;
//gotoxy(55,8);cout<<"radius  ";cin>>ObjRad;
//if (ObjRad<=0) {closegraph();return;}


int et,excitetype;
et=1;gotoxy(1,1+et);cout<<et<<". x-axis dilation";
et=2;gotoxy(1,1+et);cout<<et<<". Random 6-point 3D velocity";
et=3;gotoxy(1,1+et);cout<<et<<". Breathing";
et=4;gotoxy(1,1+et);cout<<et<<". Kick up from bottom (velocity)";

et=5;gotoxy(1,1+et);cout<<et<<". Coin cup distortion";
et=6;gotoxy(1,1+et);cout<<et<<". Coin fold distortion";
et=7;gotoxy(1,1+et);cout<<et<<". Coin saddle distortion";
et=8 ;gotoxy(1,1+et);cout<<"10+n: wet dog - order n ";
et=9 ;gotoxy(1,1+et);cout<<"20+n: Toroidal - order n ";
et=10;gotoxy(1,1+et);cout<<"30+n: Spheroidal - order n ";
et=11;gotoxy(1,1+et);cout<<"40+i: (1,2,3) axis dilation";

excitetype=44;
///gotoxy(35,6);cout<<"Excitation type? "; cin>>excitetype;
if (excitetype<0) {closegraph();return;}

et=7;gotoxy(1,1+et);cout<<"                           ";
et=8;gotoxy(1,1+et);cout<<"                           ";
et=9;gotoxy(1,1+et);cout<<"                           ";
et=10;gotoxy(1,1+et);cout<<"                           ";
et=11;gotoxy(1,1+et);cout<<"                           ";


int iClong=Clong;
int iCshear=Cshear;
//gotoxy(1,9);cout<<"Clong="<<iClong<<" m/s      ";
//gotoxy(1,10);cout<<"Cshear="<<iCshear<<" m/s    ";


float etaguess;float xiguess;
if ((excitetype>=20)&&(excitetype<30))
{ //.
gotoxy(35,7);cout<<"Eta guess? ";cin>>etaguess;
ksguess=etaguess/ObjRad;
} //.
if ((excitetype>=30)&&(excitetype<40))
{ //..
gotoxy(35,7);cout<<"Eta guess? ";cin>>etaguess;
xiguess=etaguess*Cshear/Clong;
klguess=xiguess/ObjRad;
} //..

int ticksper=4;

long int countruns=0;
int manualendflag;
float ttotperFTrun;

float binw=0.0010,binwHz,binweta,etamax,etamin;
gotoxy(35,8);cout<<"eta              ";
etamin=.4;
//gotoxy(30,8 );cout<<"eta min ";cin>>etamin;
float omeg2,r33,ftfminHz,ftfrmin,omegamin;
etamax=6.1;
//gotoxy(45,8 );cout<<"max ";cin>>etamax;
float FTpeakwidth=0.05; // width of Fourier peaks (eta)
FTpeakwidth=0.2;
//gotoxy(35,9);cout<<"peak width ";cin>>FTpeakwidth;
if (FTpeakwidth<=0) {closegraph();return;}
float MinObjRad;
//MinObjRad=3;
gotoxy(30,10);cout<<"Min R? ";cin>>MinObjRad;
if (MinObjRad<=0) {closegraph();return;}
int extrarows;

extrarows=0;
///gotoxy(42,10);cout<<"Extra rows at top? (0-8)";cin>>extrarows;
ipytop=ipytop-16*extrarows;
ftfminHz=ftfrmin/ObjRad;
int ifr,ifrcursor=nofr/2;

int ipyft;

BeginNextRun:

manualendflag=0;
countruns++;
ipyft=ipybot-countruns;

setfillstyle(SOLID_FILL,0);
bar(0,ipyft-20,639,ipyft-2);
setcolor(1);
line(0,ipytop,639,ipytop);
float eta; int ipx,ipy,ieta;

//  0.001 tick marks and labels
if (etamax-etamin<50*0.001)
for (ieta=0;ieta<1000.0*etamax;ieta=ieta+1)
{// ieta 0.01
if ((ieta%10)!=0)
{ //.
ipx = nofr*(0.001*ieta-etamin)/(etamax-etamin);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipyft-2,1);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+1,1);
if ((etamax-etamin<=0.02)&&(etamax-etamin>=0.01)&&((ieta%5)==0))
{ //..
ix=-1+1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ //...
cout.precision(2);
gotoxy(ix,iy);
if (ieta<10)
 cout<<".00"<<ieta;
else if (ieta<100)
 cout<<".0"<<ieta;
else
 cout<<0.001*ieta;
} //...
} //..
if ((etamax-etamin<=0.01)&&(etamax-etamin>=0.005)&&((ieta%2)==0))
{ //..
ix=-1+1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{
cout.precision(2);
gotoxy(ix,iy);
if (ieta<10)
 cout<<".00"<<ieta;
else if (ieta<100)
 cout<<".0"<<ieta;
else
 cout<<0.001*ieta;
}
} //..
if (etamax-etamin<0.005)
{ //..
ix=-1+1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ //...
cout.precision(2);
gotoxy(ix,iy);
if (ieta<10)
 cout<<".00"<<ieta;
else if (ieta<100)
 cout<<".0"<<ieta;
else
 cout<<0.001*ieta;
} //...
} //..
} //.
}// eta 0.001

//  0.01 tick marks and labels
if (etamax-etamin<50*0.01)
for (ieta=0;ieta<100.0*etamax;ieta=ieta+1)
{// ieta 0.01
if ((ieta%10)!=0)
{ //.
ipx = nofr*(0.01*ieta-etamin)/(etamax-etamin);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipyft-2,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+1,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+2,7);
if ((etamax-etamin<=0.2)&&(etamax-etamin>=0.1)&&((ieta%5)==0))
{ //..
ix=-1+1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ //...
cout.precision(2);
gotoxy(ix,iy);
if (ieta<10)
 cout<<".0"<<ieta;
else if (ieta<100)
 cout<<"."<<ieta;
else
 cout<<0.01*ieta;
} //...
} //..
if ((etamax-etamin<=0.1)&&(etamax-etamin>=0.05)&&((ieta%2)==0))
{ //..
ix=-1+1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{
cout.precision(2);
gotoxy(ix,iy);
if (ieta<10)
 cout<<".0"<<ieta;
else if (ieta<100)
 cout<<"."<<ieta;
else
 cout<<0.01*ieta;
}
} //..
if (etamax-etamin<0.05)
{ //..
ix=-1+1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ //...
cout.precision(2);
gotoxy(ix,iy);
if (ieta<10)
 cout<<".0"<<ieta;
else if (ieta<100)
 cout<<"."<<ieta;
else
 cout<<0.01*ieta;
} //...
} //..
} //.
}// eta 0.01

for (ieta=0;ieta<10.0*etamax;ieta=ieta+1)
{// ieta 0.1
if ((ieta%10)!=0)
{ //.
ipx = nofr*(0.1*ieta-etamin)/(etamax-etamin);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipyft-2,1);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+1,1);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+2,1);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+1,1);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+2,1);
if ((etamax-etamin<2)&&(etamax-etamin>=1)&&((ieta%5)==0))
{ //..
ix=1+(ipx-4)/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ // ...
cout.precision(1);
gotoxy(ix,iy);
if (ieta<10)
 cout<<"."<<ieta;
else
 cout<<0.1*ieta;
} // ...
} //..
if ((etamax-etamin<=1)&&(etamax-etamin>=0.5)&&((ieta%2)==0))
{ //..
ix=1+(ipx-4)/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ //
cout.precision(1);
gotoxy(ix,iy);
if (ieta<10)
 cout<<"."<<ieta;
else
 cout<<0.1*ieta;
} //
} //..
if (etamax-etamin<0.5)
{ //..
ix=1+(ipx-4)/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79-1)&&(iy>=1)&&(iy<=25))
{ //...
cout.precision(1);
gotoxy(ix,iy);
if (ieta<10)
 cout<<"."<<ieta;
else
 cout<<0.1*ieta;
} //...
} //..
} //.
}// eta 0.1
for (ieta=0;ieta<etamax;ieta=ieta+1)
{// eta
int ix,iy;
ipx = nofr*(ieta-etamin)/(etamax-etamin);
ix=1+ipx/8; iy=(ipytop+32)/16;
if ((ix>=1)&&(ix<=79)&&(iy>=1)&&(iy<=25)) {gotoxy(ix,iy);cout<<ieta;}
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipyft-2,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipyft-3,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipyft-4,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+1,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+2,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+3,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+4,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx,ipytop+5,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+1,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+2,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+3,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+4,7);
if ((ipx>=0)&&(ipx<640)) putpixel(ipx+1,ipytop+5,7);
}// ieta
//ObjRad=sqrt(14.0+countruns);
ObjRad=MinObjRad*(ipytop-1.0*ipybot)/(ipytop*1.0+countruns-1.0*ipybot);
if (ObjRad>0.98*maxpossrad)
{ // done
setrgbpalette(pal.colors[7],30,63,63); // for text
char beepc=7;
cout<<beepc;
makescr();
gotoxy(1,12);cout<<"** All done **";
getch();closegraph();
return;
} // done
cout.precision(2);
gotoxy(69,25);cout<<"R=         ";
gotoxy(69,25);cout<<"R="<<MinObjRad<<"-"<<ObjRad<<"  ";
ftfminHz=etamin*Cshear/(2.0*3.14159*ObjRad);
binw=(etamax-etamin)*Cshear/(ObjRad*nofr);
ttotperFTrun=2.0*3.14159*ObjRad/(FTpeakwidth*Cshear); // ?? check this later
//ftfminHz=ftfrmin/(2.0*3.14159*ObjRad);
//dt=dt*0.98;
//dt=0.02;
//aspectr=aspectr*1.01;
//cin>>ObjRad;
//if (ObjRad<=0) return;
//ObjRad=ObjRad+1.0;
//gotoxy(1,17);cout<<"alright at d ";
gotoxy(55,5);cout<<"shape="<<whichshape<<" ";
if (whichshape==1) {gotoxy(55,6);cout<<"radrat="<<radratio<<" ";}
if (whichshape==2) {gotoxy(55,6);cout<<"asp.r="<<aspectr<<"  ";}
if (whichshape==3) {gotoxy(55,6);cout<<"asp.r="<<aspectr<<"  ";}
gotoxy(55,8);cout<<"radius="<<ObjRad<<"  ";
//gotoxy(55,11);cout<<"dt "<<dt<<" ";
///gotoxy(55,12);cout<<"diss coef "<<dissipcoef<<" ";

mintotalenergy=2000;maxtotalenergy=-2000;
mindistever=1000;

if (drawanimation==1)
for (im=0;im<ActNm;im++)
{ // im
putpixel(bx0+oldmx[im]-128,by0-oldmy[im],0);
if (zview) putpixel(bx0+171.2+oldmz[im]-128,by0-oldmy[im],0);
} // im


for (im=0;im<MaxNm;im++)
{ // im
iatomtype[im]=0;
} // im

float angle3=0.01*random(300);
float angle13=0.01*random(300);

char crysorient[40];
//angle3=0.25*3.1415926; // for [111]
//angle13=-0.9553;       // for [111]
///strcpy(crysorient,"[111]");
angle3=0;                // for [100]
angle13=0;               // for [100]
strcpy(crysorient,"[100]");
gotoxy(30,10);cout<<"* set angle 3 "<<angle3<<" ang13="<<angle13<<" ";

ang4=0.01*random(300);
ang5=0.01*random(157);
ang4=-3.14159/2.0;
ang4=ang4+0.3; // oblique
//ang5=3.14159/2.0;
ang5=0.2;
//ang4=0;
//gotoxy(30,11);cout<<"* set angle 4 "<<ang4<<" ang5="<<ang5<<" ";

ActNm=0;

gotoxy(1,1);cout<<"Initialize positions:   ";

float extrafactor=1;
if (whichshape==0) extrafactor=1; // sphere
if (whichshape==4) extrafactor=1; // half sphere
if (whichshape==2) extrafactor=sqrt(1.0+aspectr*aspectr); //
if (whichshape==3) extrafactor=sqrt(1.0+1.0+aspectr*aspectr); //
float cosangle3=cos(angle3),sinangle3=sin(angle3);
float cosang4=cos(ang4),sinang4=sin(ang4);
float cosang5=cos(ang5),sinang5=sin(ang5);
float cosangle13=cos(angle13),sinangle13=sin(angle13);
for (iz=-1-extrafactor*ObjRad/uca;iz<=extrafactor*ObjRad/uca+1;iz++)
for (iy=-1-extrafactor*ObjRad/uca;iy<=extrafactor*ObjRad/uca+1;iy++)
for (ix=-1-extrafactor*ObjRad/uca;ix<=extrafactor*ObjRad/uca+1;ix++)
{ // im
float x1,y1,z1,x13,y13,z13,r1,r1b=ObjRad,x2,y2,z2,p866;
// Simple cubic lattice
x2=ix*uca;
y2=iy*uca;
z2=iz*uca;

x1=cosangle3*x2+sinangle3*y2;
y1=cosangle3*y2-sinangle3*x2;
z1=z2;
x13=cosangle13*x1+sinangle13*z1;
z13=cosangle13*z1-sinangle13*x1;
y13=y1;
if (whichshape==0) r1=sqrt(x13*x13+y13*y13+z13*z13); // 0. solid sphere
if (whichshape==4) r1=sqrt(x13*x13+y13*y13+z13*z13); // 4. half sphere
if (whichshape==4) r1b=sqrt(x13*x13+y13*y13+(z13-npcmz*ObjRad)*(z13-npcmz*ObjRad)); // 4. half sphere
if (whichshape==1) r1=sqrt(x13*x13+y13*y13+z13*z13);
if (whichshape==2) r1=sqrt(z13*z13+y13*y13);
int insideflag=0;
if ((whichshape==0)&&(r1<ObjRad)) insideflag=1; // solid sphere
if ((whichshape==4)&&(r1<ObjRad)&&(
 (z13<0)||(r1b<nprad*ObjRad)
                                   )) insideflag=1; // half sphere with extra part
if ((whichshape==1)&&(r1<ObjRad)&&(r1>radratio*ObjRad)) insideflag=1; // holl sphere
if ((whichshape==2)&&(r1<ObjRad)&&(fabs(x13)<aspectr*ObjRad)) insideflag=1; //cyl
if ((whichshape==3)&&(fabs(x13)<ObjRad)&&(fabs(y13)<aspectr*ObjRad)&&(fabs(z13)<ObjRad)) insideflag=1;
if (insideflag==1)
{ // in circle/cylinder/cube
if (drawanimation==1) mcol[ActNm]=14;
///if ((ix==iy)&&(iy==iz)) iatomtype[ActNm]=1;
if (r1b<=nprad*ObjRad) iatomtype[ActNm]=8;

z2=cosang4*z13+sinang4*y13;
y2=cosang4*y13-sinang4*z13;
x2=x13;
x1=cosang5*x2+sinang5*y2;
y1=cosang5*y2-sinang5*x2;
z1=z2;
mx[ActNm]=x1; my[ActNm]=y1; mz[ActNm]=z1;
if (drawanimation==1)
{ //
oldmx[ActNm]=sfx*x1+128; oldmy[ActNm]=sfx*y1;
oldmz[ActNm]=sfx*z1+128;
} //
mvx[ActNm]=ix;
mvy[ActNm]=iy;
mvz[ActNm]=iz;
ActNm++;
if (ActNm>=MaxNm)
{ // not enough memory space
gotoxy(1,1);cout<<"Error: * MaxNm = "<<MaxNm<<" is too low *"<<beepc;
//makescr();
getch();closegraph();return;
} // not enough memory space
} // in circle
} // im

int randm1,randm2,randm3;
randm1=random(ActNm);
randm2=random(ActNm);
randm3=random(ActNm);

gotoxy(1,2);cout<<"ActNm = "<<ActNm<<" "<<MaxNm<<" ";
gotoxy(72,24);cout<<"N="<<ActNm;
if (matername[0]!='?')
{ //
gotoxy(45,25);cout<<matername;
} //
else if (isotropic==1)
{ // isotropic
gotoxy(45,25);cout<<"PR="<<Poissonratio;
} // isotropic

gotoxy(1,2);cout<<"ActNm = "<<ActNm<<"  ";

gotoxy(1,3);cout<<"Find 8-clusters:";
int countnotcomplete8=0;
for (im=0;im<ActNm;im++) m9nl[im]=0;
for (im=0;im<ActNm-7;im++)
{ //.
if ((im%100)==0) {gotoxy(1,4);cout<<im<<" ";}
int ix0,iy0,iz0,ix2,iy2,iz2,count8,m11nl0;
ix0=mvx[im];
iy0=mvy[im];
iz0=mvz[im];
count8=0; m11nl0=0;
for (im2=im+1;im2<ActNm;im2++)
{ // im2
ix2=mvx[im2];
iy2=mvy[im2];
iz2=mvz[im2];
if (iz2==iz0+1)
{ // z 1
if (iy2==iy0+1)
{ // y 1
if (ix2==ix0+1) {m9nl[im]= im2-im;count8++;}
if (ix2==ix0)   {                 count8++;}
} // y 1
else if (iy2==iy0)
{ // y 0
if (ix2==ix0+1) {m11nl0   =im2-im;count8++;}
if (ix2==ix0)   {                 count8++;}
} // y 0
} // z 1
else if (iz2==iz0)
{ // z 0
if (iy2==iy0+1)
{ // y 1
if (ix2==ix0+1) {int i82=im2-im;m13nl[im]=i82;count8++;}
if (ix2==ix0)   {                 count8++;}
} // y 1
else if (iy2==iy0)
{ // y 0
if (ix2==ix0+1) {                 count8++;}
} // y 0
} // z 0
//if ((ix2==ix0+0)&&(iy2==iy0+0)&&(iz2==iz0+0)) {                 count8++;}
// m10nl = m9nl - 1
// m15nl = 1
// m16nl = 0
if (count8==8-1) break;
} // im2
m11nl0=m9nl[im]-m11nl0;
m11nl[im]=m11nl0;
if (count8!=8-1)
{ // not a complete 8-cluster
countnotcomplete8++;
m9nl[im]=0;
} // not a complete 8-cluster
} //.
gotoxy(1,4);cout<<countnotcomplete8<<" incomplete 8-clusters ";

/*
gotoxy(30,5);cout<<"Example:";
i1=random(ActNm);
if (m9nl[i1]!=0)
{ //..
int i89;
gotoxy(30,6); cout<< m9nl[i1]<<"  ";
gotoxy(30,7); cout<< m9nl[i1]-1<<"  ";       // m10 = m9 - 1
gotoxy(30,8);i89=m11nl[i1];i89=m9nl[i1]-i89;cout<<i89<<"  ";
                                                  // m11 = m9 - m11nl
gotoxy(30,9); cout<<i89-1<<"  ";                   // m12 = (m9-m11nl) - 1
gotoxy(30,10);i89=m13nl[i1];cout<<i89<<"  ";
gotoxy(30,11);i89=m13nl[i1];i89=i89-1;cout<<i89<<"  "; // m14 = m13 - 1
gotoxy(30,12);cout<<"1   ";                     // m15 = 1
gotoxy(30,13);cout<<"0   ";                     // m16 = 0
} //..
*/


gotoxy(1,3);cout<<"Making neighbour list:";

// Now construct near neighbour list:
for (im=0;im<ActNm;im++)
{ // im
m0nl[im]=0;
m1nl[im]=0;
m2nl[im]=0;
m3nl[im]=0;
m4nl[im]=0;
m5nl[im]=0;
m6nl[im]=0;
m7nl[im]=0;
m8nl[im]=0;
} //
int maxnoneigh=0;
for (im1=0;im1<ActNm;im1++)
{ // im1
int mnn1=0,mnn2=6/2;//,mnn3=(6+12)/2;
if ((im1%100)==0) {gotoxy(1,4);cout<<im1<<" ";}
for (im2=0;im2<im1;im2++)
{ // im2
dx=(mx[im1]-mx[im2]);
dy=(my[im1]-my[im2]);
dz=(mz[im1]-mz[im2]);
float d2;
d2=dx*dx+dy*dy+dz*dz;
//if (d2<1.5*1.5) // for nearest neighbours only
//if (d2<2.6*2.6)
if ((d2>uca*uca*0.9*0.9)&&(d2<uca*uca*1.1*1.1))
{ //..1
if (mnn1<MaxNn)
{ //.
setmnl(im1,mnn1,im1-im2);
mnn1++;
if (mnn1>maxnoneigh) {maxnoneigh=mnn1;}
if (mnn1>3) {cout<<"mnn1 reached 3 ";exit(0);}
} //.
} //..1
if ((d2>uca*uca*1.3*1.3)&&(d2<uca*uca*1.5*1.5))
{ //..1
if (mnn2<MaxNn)
{ //.
setmnl(im1,mnn2,im1-im2);
mnn2++;
if (mnn2>maxnoneigh) {maxnoneigh=mnn2;}
if (mnn2>(6+12)/2) {cout<<"mnn2 reached 9 ";exit(0);}
} //.
} //..1
} // im2
} // im1

///gotoxy(1,8);cout<<"Max no neigh = "<<maxnoneigh<<" ";

gotoxy(1,1);cout<<"Starting elast... with "<<ActNm<<" "<<MaxNm<<" ";

ttot=0;

long int litcount=0,litcountlast=0;

if (ActNm>=MaxNm)
{ //..
cout<<"*** Error ActNm > MaxNm ***";getch();return;
} //..

long int liclocklast,liclock2,liclocklastdrawn;
liclocklast=clock();
liclocklastdrawn=clock();
int graphflag=1;

for (ifr=0;ifr<nofr;ifr++)
{ // ifr
fsumr1[ifr]=0;fsumi1[ifr]=0;
fsumr2[ifr]=0;fsumi2[ifr]=0;
fsumr3[ifr]=0;fsumi3[ifr]=0;
oldfgraf[ifr]=0;
} // ifr

long int lipx,lipy,atomtowatch1,atomtowatch2,atomtowatch3;
LB4001: atomtowatch1=random(ActNm);
if ((iatomtype[atomtowatch1]!=8)&&(!kbhit())) goto LB4001;
LB4002: atomtowatch2=random(ActNm);
if ((iatomtype[atomtowatch2]!=8)&&(!kbhit())) goto LB4002;
LB4003: atomtowatch3=random(ActNm);
if ((iatomtype[atomtowatch3]!=8)&&(!kbhit())) goto LB4003;

// Set all initial velocities to zero
for (im=0;im<ActNm;im++)
{ // im
mvx[im]=0;
mvy[im]=0;
mvz[im]=0;
//iatomtype[im]=0;
} // im

//  Jan 23 2003
for (im=0;im<ActNm;im++)
{ // im Jan 23 2003
///if (mz[im]<-0.5*ObjRad) iatomtype[im]=7;
float r2;
r2=mx[im]*mx[im]+my[im]*my[im]+mz[im]*mz[im];
if (iatomtype[im]!=8)
{ //.
if (r2>(ObjRad-1)*(ObjRad-1)) iatomtype[im]=6; // 6=damped atom
} //.
} // im Jan 23 2003
//gotoxy(1,ipytop/16);cout<<VERN<<"  "<<crysorient<<
//" oriented nanosphere r="<<nprad<<"R, cmz="<<npcmz<<"R          scale is eta";

if (whichshape==0)
{ //..
gotoxy(28,25);cout<<"Sphere";
} //..
if (whichshape==2)
{ //
gotoxy(28,25);cout<<"Cyl.";
gotoxy(33,25);cout<<"h/D";
gotoxy(37,25);cout<<"=";
gotoxy(39,25);cout<<aspectr;
} //
if (whichshape==4)
{ //..
gotoxy(28,24);cout<<"Sphere on";
gotoxy(28,25);cout<<"hemisphere";
} //..

if (excitetype==1)
{ // 1
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"x-axis dilation";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]!=7) mx[im]=mx[im]*1.001;
} // im
} // 1

if (excitetype==41)
{ // 41
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"x-axis dilation";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]!=7) mx[im]=mx[im]*1.001;
} // im
} // 41

if (excitetype==44)
{ // 41
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"np x-axis dilat. !!! exag";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]==8) mx[im]=mx[im]*1.07;
} // im
} // 44

if (excitetype==46)
{ // 46
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"nsph z-axis dilat.";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]==8) mz[im]=mz[im]*1.001;
} // im
} // 46

if (excitetype==47)
{ // 47
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"np (x+y)-axis dilat.";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]==8) mx[im]=mx[im]*1.001;
if (iatomtype[im]==8) my[im]=my[im]*1.001;
} // im
} // 47

if (excitetype==42)
{ // 42
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"y-axis dilation";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]!=7) my[im]=my[im]*1.001;
} // im
} // 42

if (excitetype==43)
{ // 43
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"z-axis dilation";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]!=7) mz[im]=mz[im]*1.001;
} // im
} // 43

if (excitetype==2)
{ // 2
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"Random 6 pt.";
int ir5;
ir5=random(ActNm);mvx[ir5]-=0.01*Cshear;
ir5=random(ActNm);mvx[ir5]+=0.01*Cshear;
ir5=random(ActNm);mvy[ir5]-=0.01*Cshear;
ir5=random(ActNm);mvy[ir5]+=0.01*Cshear;
ir5=random(ActNm);mvz[ir5]-=0.01*Cshear;
ir5=random(ActNm);mvz[ir5]+=0.01*Cshear;
} // 2

if (excitetype==201)
{ // 2
gotoxy(2,23);cout<<"Init.";gotoxy(8,23);cout<<"Excit.";
gotoxy(2,24);cout<<"Random 6 pt.";
gotoxy(2,25);cout<<"in nanosphere";
int ir5;
LB201a: ir5=random(ActNm); if (iatomtype[ir5]!=8) goto LB201a;
mvx[ir5]-=0.01*Cshear;
LB201b: ir5=random(ActNm); if (iatomtype[ir5]!=8) goto LB201b;
mvx[ir5]+=0.01*Cshear;
LB201c: ir5=random(ActNm); if (iatomtype[ir5]!=8) goto LB201c;
mvy[ir5]-=0.01*Cshear;
LB201d: ir5=random(ActNm); if (iatomtype[ir5]!=8) goto LB201d;
mvy[ir5]+=0.01*Cshear;
LB201e: ir5=random(ActNm); if (iatomtype[ir5]!=8) goto LB201e;
mvz[ir5]-=0.01*Cshear;
LB201f: ir5=random(ActNm); if (iatomtype[ir5]!=8) goto LB201f;
mvz[ir5]+=0.01*Cshear;
} // 2

if (excitetype==3)
{ // 3
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"Breathing";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]!=7) mx[im]=mx[im]*1.001;
if (iatomtype[im]!=7) my[im]=my[im]*1.001;
if (iatomtype[im]!=7) mz[im]=mz[im]*1.001;
} // im
} // 3

if (excitetype==4)
{ // 4
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"Bottom cap kick";
for (im=0;im<ActNm;im++)
{ // im
if (my[im]<-(ObjRad-1.5)) mvy[im]=0.01;
} // im
} // 4

if (excitetype==7)
{ //
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"Cup";
for (im=0;im<ActNm;im++)
{ // im
mx[im]=mx[im]+0.001*(my[im]*my[im]+mz[im]*mz[im]);
} // im
} //

if (excitetype==8)
{ // 8
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"zx Folding";
for (im=0;im<ActNm;im++)
{ // im
mx[im]=mx[im]+0.001*mz[im];
} // im
} // 8

if (excitetype==801)
{ // 801
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"nsph zx Folding";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]==8) mx[im]=mx[im]+0.001*mz[im];
} // im
} // 801

if (excitetype==9)
{ // 8
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"Saddle";
for (im=0;im<ActNm;im++)
{ // im
mx[im]=mx[im]+0.001*(my[im]*my[im]-mz[im]*mz[im]);
} // im
} // 8

if ((excitetype>=10)&&(excitetype<20))
{ // my naive toroidal mode of order n
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
if (excitetype!=11) {gotoxy(2,25);cout<<"n="<<excitetype-10<<" naive toroidal";}
if (excitetype==11) {gotoxy(2,25);cout<<"Wet dog";}
for (im=0;im<ActNm;im++)
{ // im
float costheta,r,legendre;
r=sqrt(mx[im]*mx[im]+my[im]*my[im]+mz[im]*mz[im]);
if (r!=0)
{ //..
costheta=mz[im]/r;
float r5=LegenP(excitetype-10,costheta);
mvx[im]= 0.00001*Cshear*r5*r*my[im];  // July 2 added factor of r
mvy[im]=-0.00001*Cshear*r5*r*mx[im];  // to get original "wet dog" mode back
mvz[im]=0;  // Jul 30 added factor of Cshear
} //..
else
{ // origin case
mvx[im]=0;
mvy[im]=0;
mvz[im]=0;
} // origin case
} // im
} // 6

if (excitetype==111)
{ // z twist
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"nsph z twist";
for (im=0;im<ActNm;im++)
{ // im
if (iatomtype[im]==8)
{ // ==8
mvx[im]= 0.00001*Cshear*mz[im]*my[im];
mvy[im]=-0.00001*Cshear*mz[im]*mx[im];
mvz[im]=0;
} // ==8
} // im
} // z twist

if ((excitetype>=20)&&(excitetype<30))
{ // true toroidal mode of order n
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"n="<<excitetype-20<<" toroidal";
for (im=0;im<ActNm;im++)
{ // im
mvx[im]=0.00001*Cshear*CrPix(excitetype-20,mx[im],my[im],mz[im]);
mvy[im]=0.00001*Cshear*CrPiy(excitetype-20,mx[im],my[im],mz[im]);
mvz[im]=0.00001*Cshear*CrPiz(excitetype-20,mx[im],my[im],mz[im]);
//  Jul 31  added factor of Cshear; changed C2rPi to CrPi
} // im
} // true toroidal mode of order n

if ((excitetype>=30)&&(excitetype<40))
{ // true spheroidal mode of order n
gotoxy(2,24);cout<<"Init.";gotoxy(8,24);cout<<"Excit.";
gotoxy(2,25);cout<<"n="<<excitetype-30<<" spheroidal";
for (im=0;im<ActNm;im++)
{ // im
mvx[im]=0.0001*Cshear*THETAx(excitetype-30,mx[im],my[im],mz[im]);
mvy[im]=0.0001*Cshear*THETAy(excitetype-30,mx[im],my[im],mz[im]);
mvz[im]=0.0001*Cshear*THETAz(excitetype-30,mx[im],my[im],mz[im]);
//  Jul 31  added factor of Cshear
} // im
} // true spheroidal mode of order n

float KEatstart=0;
for (im=0;im<ActNm;im++)
{ // im
float v2now;
cellmass=densityrho*uca*uca*uca;
v2now=(mvx[im]*mvx[im]+mvy[im]*mvy[im]+mvz[im]*mvz[im]);
KEatstart=KEatstart+0.5*cellmass*v2now;
} // im

makescr(); // initial blank frame


LBnexttimestep:


if (ttot>ttotperFTrun)
{ //
///goto BeginNextRun;
fclose(scrfout);
return;
} //

if ((litcount%5)==0)
{ //....
//sfx=sfx*1.2;
setfillstyle(SOLID_FILL,0);
int bx0=100;
int by0=200;
bar(bx0-85-1,by0-85-1,bx0+85+1,by0+85+1);
setcolor(7);
rectangle(bx0-85-1,by0-85-1,bx0+85+1,by0+85+1);
float zlayer;
for (zlayer=-ObjRad;zlayer<ObjRad+0.2;zlayer=zlayer+0.2)
{ // zlayer
for (im=0;im<ActNm;im++)
{ // im
int ic;
int i54=sfx*mx[im]+128;
int i55=sfx*my[im],ic5=mcol[im];
int i56=sfx*mz[im]+128;
ic=mcol[im];
//if ((ic==15)&&(mz[im]>zlayer-1.0)&&(mz[im]<=zlayer))
if ((mcol[im]!=0)&&(mz[im]>zlayer-1.0)&&(mz[im]<=zlayer))
{ //
int ic2,iring,noring=14;
for (iring=0;iring<noring;iring++)
{ // iring
ic2=ic-iring;if (ic2<1) ic2=1;
if (ic2==7) ic2=6; // blue
if (ic2==1) ic2=2; // red
circleb(bx0+sfx*mx[im],by0-sfx*my[im],iring*sfx*0.5*(1.4142136*uca)/noring,(iring+1)*sfx*0.5*(1.4142136*uca)/noring,ic2);
} // iring
} //
} // im
if (kbhit()) {getch(); return;}
} // zlayer
makescr();
licountframes++;
gotoxy(40,16);cout<<licountframes<<" frames ";
//closegraph();
//return;
} //....

litcount++;
ttot=ttot+dt;

for (ifr=0;ifr<nofr;ifr++)
{ // .. ifr
omeg2=(2.0*3.1415926*ftfminHz)+ifr*binw;
float rr,ri;
rr=cos(omeg2*ttot);ri=sin(omeg2*ttot);
fsumr1[ifr]+=mvx[atomtowatch1]*rr;
fsumi1[ifr]+=mvx[atomtowatch1]*ri;
fsumr2[ifr]+=mvy[atomtowatch2]*rr;
fsumi2[ifr]+=mvy[atomtowatch2]*ri;
fsumr3[ifr]+=mvz[atomtowatch3]*rr;
fsumi3[ifr]+=mvz[atomtowatch3]*ri;
} // .. ifr

if (dissipcoef!=0)
{ //..
for (im1=0;im1<ActNm;im1++)
{ // im1
max[im1]=0;may[im1]=0;maz[im1]=0;
} // im1
} //..

float eps=0.0174989; // 6-12 energy parameter (makes dF/dr=ksp=1.0)

float mmass=1.0; // mass of a molecule
float d1x,d1y,d1z,d2,dp1x,dp1y,dp1z;

float countdamagedatoms=0;

for (im1=0;im1<ActNm;im1++)
{ // im1
d1x=mx[im1]; d1y=my[im1]; d1z=mz[im1];
dp1x=0; dp1y=0; dp1z=0;
int i5,dim12,i1n;
float maxr2fordamage=0;

if (m9nl[im1]!=0)
{ // if atom im1 is corner of a cluster
int i8,i81; float cm8x,cm8y,cm8z,sum8,dx,dy,dz,dpx,dpy,dpz;
n8pt[0]=im1+m9nl[im1];
n8pt[1]=n8pt[0]-1;
i81=m11nl[im1];
i81=m9nl[im1]-i81;
n8pt[2]=im1+i81;
n8pt[3]=n8pt[2]-1;
i81=m13nl[im1];
n8pt[4]=im1+i81;
n8pt[5]=n8pt[4]-1;
n8pt[6]=im1+1;
n8pt[7]=im1;
cm8x=0;
cm8y=0;
cm8z=0;
for (i8=0;i8<8;i8++)
{ // i8
cm8x+=mx[n8pt[i8]];
cm8y+=my[n8pt[i8]];
cm8z+=mz[n8pt[i8]];
} // i8
cm8x=cm8x/8.0;
cm8y=cm8y/8.0;
cm8z=cm8z/8.0;
sum8=0;
for (i8=0;i8<8;i8++)
{ // i8
dx=mx[n8pt[i8]]-cm8x;
dy=my[n8pt[i8]]-cm8y;
dz=mz[n8pt[i8]]-cm8z;
sum8=sum8+dx*dx+dy*dy+dz*dz;
} // i8
sum8=sum8-6.0*uca*uca;
sum8=sum8*k8pt/(12.0*uca*uca);
cellmass=densityrho*uca*uca*uca;
sum8=sum8*dt/cellmass;
for (i8=0;i8<8;i8++)
{ // i8
mvx[n8pt[i8]]+=(cm8x-mx[n8pt[i8]])*sum8;
mvy[n8pt[i8]]+=(cm8y-my[n8pt[i8]])*sum8;
mvz[n8pt[i8]]+=(cm8z-mz[n8pt[i8]])*sum8;
} // i8
} // if atom im1 is corner of a cluster

// ksp1 and ksp2 2-body forces:
for (i1n=0;i1n<MaxNn;i1n++)
{ // i1
float Fbyr,Fx,Fy,Fz,r121;
dim12=getmnl(im1,i1n);
if (dim12!=0)
{ // dim12!=0
im2=im1-dim12;
dx=(d1x-mx[im2]); dy=(d1y-my[im2]); dz=(d1z-mz[im2]);
r122=dx*dx+dy*dy+dz*dz;
//if (r122>maxr2fordamage) maxr2fordamage=r122;
//if (r122>maxeverr122) maxeverr122=r122;
//if (r122<mineverr122) mineverr122=r122;
//if (r122>maxr122now) maxr122now=r122;
//if (r122<minr122now) minr122now=r122;
if (r122>0.01*uca*uca)
{ // not too close
float vr,Fvbyr,Fvx,Fvy,Fvz;
if (LennJonesflag)
{ //.. // Lennard-Jones 6-12 potential
// U(r) = 4.0*eps*(1/r^12-1/r^6)
// F(r) = dU/dr = 4.0*eps*(12/r^13-6/r^7)
// F(r)/r = 4.0*eps*(12/r^14-6/r^8...
// dF/dr = 4.0*eps*(156/r^14-42/r^8)
//         4.0*eps*(30.9543205-16.667711)=57.14644*eps = 1.0
r128=r122*r122*r122*r122;
r1214=r128*r122*r122*r122;
Fbyr=4.0*eps*((12.0/r1214)-(6.0/r128));
} //..
else
{ //.. // Quadratic potential;
r121=sqrt(r122);
if (i1n<6/2)
Fbyr=-ksp1*(r121-uca)/r121;
if ((i1n>=6/2)&&(i1n<(6+12)/2))
Fbyr=-ksp2*(r121-1.4142136*uca)/r121;
} //..
Fx=dx*Fbyr; Fy=dy*Fbyr; Fz=dz*Fbyr;
dp1x+=Fx; dp1y+=Fy; dp1z+=Fz;

if (dissipcoef==0)
{ //.
if (iatomtype[im2]!=7)
{ // Jan 23 2003
cellmass=densityrho*uca*uca*uca;
mvx[im2]-=dt*Fx/cellmass;
mvy[im2]-=dt*Fy/cellmass;
mvz[im2]-=dt*Fz/cellmass;
} // Jan 23 2003
} //.
else
{ //..
cellmass=densityrho*uca*uca*uca;
max[im2]-=Fx/cellmass; may[im2]-=Fy/cellmass; maz[im2]-=Fz/cellmass;
vr=dx*(mvx[im1]-mvx[im2])+dy*(mvy[im1]-mvy[im2]);
vr=vr;
Fvbyr=-dissipcoef*vr/r128;
Fvx=Fvbyr*dx; Fvy=Fvbyr*dy; Fvz=Fvbyr*dz;
max[im1]+=Fvx/cellmass; may[im1]+=Fvy/cellmass; maz[im1]+=Fvz/cellmass;
max[im2]-=Fvx/cellmass; may[im2]-=Fvy/cellmass; maz[im2]-=Fvz/cellmass;
} //..

} // not too close
} // dim12!=0
} // i1
if (dissipcoef==0)
{ //..
if (iatomtype[im1]!=7)
{ // Jan 23 2003
cellmass=densityrho*uca*uca*uca;
mvx[im1]+=dt*dp1x/cellmass;
mvy[im1]+=dt*dp1y/cellmass;
mvz[im1]+=dt*dp1z/cellmass;
} // Jan 23 2003
} //..
else
{ //.
cellmass=densityrho*uca*uca*uca;
max[im1]+=dp1x/cellmass;
may[im1]+=dp1y/cellmass;
maz[im1]+=dp1z/cellmass;
} //.
if (maxr2fordamage>3.5) countdamagedatoms=countdamagedatoms+1.0;
} // im1

if (dissipcoef!=0)
{ //.
for (im=0;im<ActNm;im++)
{ // ix iy
mvx[im]=mvx[im]+dt*max[im];
mvy[im]=mvy[im]+dt*may[im];
mvz[im]=mvz[im]+dt*maz[im];
} // ix iy
} //.

{ // Jan 25, 2003 simple surface damping
// This is a simple kind of surface damping which is matched to
// long wavelength longitudinal plane waves.  It is not adjusted
// to efficiently absorb transverse plane waves.
// Rough experimental optimization of dampadjustfac for random
// point disturbances in a sphere is 0.5.  cc3bdmp1.cpp
float dampfac,Clong,cellmass,bdamp;
Clong=sqrt(C11/densityrho);
cellmass=(densityrho*uca*uca*uca);
bdamp=densityrho*Clong*uca*uca;
dampfac=(1.0-dampadjustfac*dt*bdamp/cellmass);
// Surface damping feature:
for (im=0;im<ActNm;im++)
{ //.
if (iatomtype[im]==6)
{ //..
mvx[im]=mvx[im]*dampfac;
mvy[im]=mvy[im]*dampfac;
mvz[im]=mvz[im]*dampfac;
} //..
} //.
} //

for (im=0;im<ActNm;im++)
{ // ix iy
if (iatomtype[im]!=7)
{ // Jan 23 2003
mx[im]=mx[im]+dt*mvx[im];
my[im]=my[im]+dt*mvy[im];
mz[im]=mz[im]+dt*mvz[im];
} // Jan 23 2003
} // ix iy

if (clock()>=liclocklastdrawn+ticksper)
{ //. another clock tick

gotoxy(39,1);cout<<litcount<<" ";
long int microsecperit=ticksper*1.0e6/((litcount-litcountlast)*CLK_TCK);
gotoxy(39,3);cout<<" "<<microsecperit<<" æs/it ";

float KEnow=0,peakv2now=0,v2now;
for (im=0;im<ActNm;im++)
{ // im
v2now=(mvx[im]*mvx[im]+mvy[im]*mvy[im]+mvz[im]*mvz[im]);
KEnow=KEnow+0.5*cellmass*v2now;
if (v2now>peakv2now) peakv2now=v2now;
} // im
cout.precision(2);
gotoxy(1,5);cout<<"KE="<<KEnow<<"   ";
if (KEatstart!=0)
{ //...
cout.precision(3);
gotoxy(20,5);cout<<" KE/KE = "<<KEnow/KEatstart<<"   ";
//if (0.5*litcount<640)
//putpixel(0.5*litcount,480-540*KEnow/KEatstart,dampadjustcol);
} //...
gotoxy(1,6);cout<<"KE/atom = "<<KEnow/ActNm<<"  ";
//gotoxy(1,7);cout<<"rms v = "<<sqrt(2.0*(KEnow/ActNm))<<"  ";
gotoxy(1,7);cout<<"peak v = "<<sqrt(peakv2now)<<"  ";


liclocklastdrawn=clock();
litcountlast=litcount;

if (drawanimation==1)
for (im=0;im<ActNm;im++)
{ // ix iy
if (mcol[im]!=0)
{ //
putpixel(bx0+oldmx[im]-128,by0-oldmy[im],0);
if (zview) putpixel(bx0+171.2+oldmz[im]-128,by0-oldmy[im],0);
}
} // ix iy

for (im=0;im<ActNm;im++)
{ // ix iy
// This part is to avoid program crashing in graphics calls
if (mx[im]<-100) mx[im]=-100;
if (my[im]<-100) my[im]=-100;
if (mz[im]<-100) mz[im]=-100;
if (mx[im]>100) mx[im]=100;
if (my[im]>100) my[im]=100;
if (mz[im]>100) mz[im]=100;

if (drawanimation==1)
if (mcol[im]!=0)
{ // non-black
int ic; ic=14;
int i54=sfx*mx[im]+128;
int i55=sfx*my[im];
int i56=sfx*mz[im]+128;
float r1,r2,dx,dy,dz;

ic=15; //
if (iatomtype[im]==7) ic=7; // Jan 23 2003  frozen atom
if (iatomtype[im]==6) ic=1; // Jan 25 2003  damped atom
if (iatomtype[im]==8) ic=7; // Jan 25 2003  nanoparticle atom

if (drawanimation)
if (ic!=0)
{ //
putpixel(bx0+i54-128,by0-i55,ic);
if (zview) putpixel(bx0+171.2+i56-128,by0-i55,ic);
///if (iatomtype[im]==1) {setcolor(1);circle(bx0+i54-128,by0-i55,2);}
///if (iatomtype[im]==1) {setcolor(1);circle(bx0+171.2+i56-128,by0-i55,2);}
} //
} // non-black
if (drawanimation==1)
{ //
oldmx[im]=mx[im]*sfx+128;
oldmy[im]=my[im]*sfx;
oldmz[im]=mz[im]*sfx+128;
} //
} // ix iy

float fhz;
/*
for (fhz=0.001;fhz<0.5;fhz=fhz+0.001)
{ // fhz
omeg2=2.0*3.14159*ftfminHz+fhz*2.0*3.14159;
i1=omeg2/binw;
putpixel(i1,200,10);
putpixel(i1,479,10);
} // fhz
*/
/*
for (fhz=0.01;fhz<0.5;fhz=fhz+0.01)
{ // fhz
omeg2=2.0*3.14159*ftfminHz+fhz*2.0*3.14159;
i1=omeg2/binw;
putpixel(i1,200,13);
putpixel(i1,479,13);
} // fhz
for (fhz=0.1;fhz<0.5;fhz=fhz+0.1)
{ // fhz
omeg2=2.0*3.14159*ftfminHz+fhz*2.0*3.14159;
i1=omeg2/binw;
putpixel(i1,479,14);
putpixel(i1,478,14);
putpixel(i1,477,14);
} // fhz
*/

float maxr332=0;
int i8,i80;
/*
setcolor(0);
i80=0;
for (ifr=1;ifr<nofr;ifr++)
{ // .. ifr
i8=oldfgraf[ifr];
line(ifr-1,477-i80,ifr,477-i8);
i80=i8;
} // .. ifr
*/

i1=0; if (etamin==0) i1=6;
for (ifr=i1;ifr<nofr;ifr++)
{ // .. ifr
omeg2=(2.0*3.14159*ftfminHz)+ifr*binw;
r33=fsumr1[ifr]*fsumr1[ifr]+fsumi1[ifr]*fsumi1[ifr]+
    fsumr2[ifr]*fsumr2[ifr]+fsumi2[ifr]*fsumi2[ifr]+
    fsumr3[ifr]*fsumr3[ifr]+fsumi3[ifr]*fsumi3[ifr];
if (r33>maxr332) {maxr332=r33;}
} // .. ifr
if (maxr332<=0) maxr332=1.0;

//setcolor(1);line(ifrcursor,471,ifrcursor,479);
//cout.precision(4);
//gotoxy(1,10);cout<<(2.0*3.14159/Cshear)*(ifrcursor*binw/(2.0*3.1416)+ftfminHz)*ObjRad<<" ";

/*
i80=0;
for (ifr=1;ifr<nofr;ifr++)
{ // ..
r33=fsumr1[ifr]*fsumr1[ifr]+fsumi1[ifr]*fsumi1[ifr]+
    fsumr2[ifr]*fsumr2[ifr]+fsumi2[ifr]*fsumi2[ifr]+
    fsumr3[ifr]*fsumr3[ifr]+fsumi3[ifr]*fsumi3[ifr];
r33=sqrt(fabs(r33)/maxr332);
setcolor(11);
i8=221.0*r33;
line(ifr-1,477-i80,ifr,477-i8);
oldfgraf[ifr]=i8;
i80=i8;
} // ..
*/


for (ifr=0;ifr<nofr;ifr++)
{ // ..
r33=fsumr1[ifr]*fsumr1[ifr]+fsumi1[ifr]*fsumi1[ifr]+
    fsumr2[ifr]*fsumr2[ifr]+fsumi2[ifr]*fsumi2[ifr]+
    fsumr3[ifr]*fsumr3[ifr]+fsumi3[ifr]*fsumi3[ifr];
r33=sqrt(fabs(r33)/maxr332);
ic=15*r33; if (ic>15) ic=15;
if (ic==7) ic=8;
if (ic==1) ic=2;
putpixel(ifr,ipyft,ic);
} // ..

LBcheckforkeystroke:
if (kbhit())
{ // key hit
char ch1;
ch1=getch();
if (ch1=='A') {dampadjustfac=dampadjustfac*1.25;gotoxy(40,15);cout<<" "<<dampadjustfac<<" ";dampadjustcol=8-dampadjustcol;}
if (ch1=='a') {dampadjustfac=dampadjustfac*0.8;gotoxy(40,15);cout<<" "<<dampadjustfac<<" ";dampadjustcol=8-dampadjustcol;}
if (ch1=='q') goto QUIT;
if (ch1=='r') makescr();
if (ch1=='a') goto BeginNextRun;
if (ch1=='t') {ticksper=ticksper+1;gotoxy(40,5);cout<<ticksper<<" ";}
if (ch1=='T') {ticksper=ticksper-1;gotoxy(40,5);cout<<ticksper<<" ";}
if (ch1=='f') {dt=dt*1.2;gotoxy(40,4);cout<<" dt="<<dt<<"  ";}
if (ch1=='e') manualendflag=1;
if (ch1=='s') {dt=dt/1.2;gotoxy(40,4);cout<<" dt="<<dt<<"  ";}

if (ch1=='c') cleardevice();
goto LBcheckforkeystroke;
} // key hit

cmx=0; cmy=0; cmz=0;
cmvx=0;cmvy=0;cmvz=0;
float ybottom=10000;

for (im=0;im<ActNm;im++)
{ // im
cmx=cmx+mx[im];cmy=cmy+my[im];cmz=cmz+mz[im];
cmvx=cmvx+mvx[im];cmvy=cmvy+mvy[im];cmvz=cmvz+mvz[im];
if (my[im]<ybottom) {ybottom=my[im];}
} // im
cmx=cmx/ActNm;cmy=cmy/ActNm;cmz=cmz/ActNm;
cmvx=cmvx/ActNm;cmvy=cmvy/ActNm;cmvz=cmvy/ActNm;
cout.precision(7);

} //.  another clock tick


goto LBnexttimestep;

QUIT:
closegraph();
fclose(scrfout);
return;
} // main

void setmnl(int im1,int in2,int im3)
{ // setmnl                                                setmnl
if (in2>=MaxNn) {cout<<"* Error: Call to setmnl() has in2="<<in2<<" ";exit(0);}
if (in2==0) {m0nl[im1]=im3; return;}
if (in2==1) {m1nl[im1]=im3; return;}
if (in2==2) {m2nl[im1]=im3; return;}
if (in2==3) {m3nl[im1]=im3; return;}
if (in2==4) {m4nl[im1]=im3; return;}
if (in2==5) {m5nl[im1]=im3; return;}
if (in2==6) {m6nl[im1]=im3; return;}
if (in2==7) {m7nl[im1]=im3; return;}
if (in2==8) {m8nl[im1]=im3; return;}
cout<<"***Error in setmnl***";exit(0);
} // setmnl                                                setmnl

int getmnl(int im1,int in2)
{ // getmnl                                                getmnl
if (in2>=MaxNn) {cout<<"* Error: Call to getmnl() has in2="<<in2<<" ";exit(0);}
int im3;
if (in2==0) {im3=m0nl[im1];return im3;}
if (in2==1) {im3=m1nl[im1];return im3;}
if (in2==2) {im3=m2nl[im1];return im3;}
if (in2==3) {im3=m3nl[im1];return im3;}
if (in2==4) {im3=m4nl[im1];return im3;}
if (in2==5) {im3=m5nl[im1];return im3;}
if (in2==6) {im3=m6nl[im1];return im3;}
if (in2==7) {im3=m7nl[im1];return im3;}
if (in2==8) {im3=m8nl[im1];return im3;}
cout<<"***Error in getmnl***";exit(0);
} // getmnl                                                getmnl

float LegenP(int n,float x)
{ //  LegenP - Legendre polynomials
float p;
if (n==0) return 1;
if (n==1) return x;
if (n==2) return 0.5*(3.0*x*x-1.0);
if (n==3) return 0.5*(5.0*x*x*x-3.0*x);
if (n==4) return 0.125*(35.0*x*x*x*x-30.0*x*x+3);
if (n==5) return 0.125*(63.0*x*x*x*x*x-70.0*x*x*x+15*x);
cout<<"* Legendre polynomial order not available *";exit(0);
} //  LegenP

float SphBess1j(int n,float x)
{ //  SphBess1j
if (n==0) return sin(x)/x;
if (n==1) return (sin(x)/(x*x))-(cos(x)/x);
if (n==2) return ( (3.0/(x*x*x)) - (1.0/x) )*sin(x) - (3.0/(x*x))*cos(x);
return SeriesSphBess1j(n,x); // use series expression (valid x<14)
///cout<<"* Spherical Bessel function order not available *";exit(0);
} //  SphBess1j


float rPix(int n,float x,float y,float z)
{ //  rPi x-component
float r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{
costheta=z/r;
return x*SphBess1j(n,ksguess*r)*LegenP(n,costheta);
}
else return 0;
} //  r Pi x-component

float rPiy(int n,float x,float y,float z)
{ //  r Pi y-component
float r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
costheta=z/r;
return y*SphBess1j(n,ksguess*r)*LegenP(n,costheta);
} //
else return 0;
} //  r Pi y-component

float rPiz(int n,float x,float y,float z)
{ //  r Pi z-component
float r,costheta;
r=sqrt(x*x+y*y+z*z);
if (r!=0)
{ //
costheta=z/r;
return z*SphBess1j(n,ksguess*r)*LegenP(n,costheta);
} //
else return 0;
} //  r Picz-component

float CrPix(int n,float x,float y,float z)
{ //  CurlrPix
float x1,x2,y1,y2,z1,z2,dx=0.01,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (rPiy(n,x,y,z2)-rPiy(n,x,y,z1))/dx;
ans -= (rPiz(n,x,y2,z)-rPiz(n,x,y1,z))/dx;
return ans;
} //  CurlrPix

float CrPiy(int n,float x,float y,float z)
{ //  CurlrPiy
float x1,x2,y1,y2,z1,z2,dx=0.01,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (rPiz(n,x2,y,z)-rPiz(n,x1,y,z))/dx;
ans -= (rPix(n,x,y,z2)-rPix(n,x,y,z1))/dx;
return ans;
} //  CurlrPiy

float CrPiz(int n,float x,float y,float z)
{ //  CrPiz
float x1,x2,y1,y2,z1,z2,dx=0.01,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (rPix(n,x,y2,z)-rPix(n,x,y1,z))/dx;
ans -= (rPiy(n,x2,y,z)-rPiy(n,x1,y,z))/dx;
return ans;
} //  CrPiz

//--------------------

float C2rPix(int n,float x,float y,float z)
{ //  curl of curl of rPi - x component
float x1,x2,y1,y2,z1,z2,dx=0.01,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (CrPiy(n,x,y,z2)-CrPiy(n,x,y,z1))/dx;
ans -= (CrPiz(n,x,y2,z)-CrPiz(n,x,y1,z))/dx;
return ans;
} //  curl of curl of rPi - x component

float C2rPiy(int n,float x,float y,float z)
{ //  CurlrPiy
float x1,x2,y1,y2,z1,z2,dx=0.01,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (CrPiz(n,x2,y,z)-CrPiz(n,x1,y,z))/dx;
ans -= (CrPix(n,x,y,z2)-CrPix(n,x,y,z1))/dx;
return ans;
} //  CurlrPiy

float C2rPiz(int n,float x,float y,float z)
{ //  C2rPiz
float x1,x2,y1,y2,z1,z2,dx=0.01,ans;
x1=x-0.5*dx;x2=x+0.5*dx;
y1=y-0.5*dx;y2=y+0.5*dx;
z1=z-0.5*dx;z2=z+0.5*dx;
ans =  (CrPix(n,x,y2,z)-CrPix(n,x,y1,z))/dx;
ans -= (CrPiy(n,x2,y,z)-CrPiy(n,x1,y,z))/dx;
return ans;
} //  C2rPiz

double factorial(double arg)
{ //
double ans,arg2,r2; int i1,i2;
ans=1.0;
i1 = floor(arg+0.5);
for (i2=2;i2<=i1;i2++)
{ // i2
r2=i2;
ans=ans*r2;
} // i2
return ans;
} //

float SeriesSphBess1j(int n,float x)
{ // infinite series spherical bessel function of 1st kind
// Note: summing to 15 gives optimal precision out to
//  about x=15, after which the series blow up, if "float"
//  is used.  If float precision is used, then it is possible
//  to sum more terms to extend the range of accurate convergence
//  without getting the blow up problem.
// (Double precision as used here is not really necessary.)
if (x>14) {cout<<"x ="<<x<<" too large for series sph bess func.";exit(0);}
float ans;
float ansd=0,nr,sr,xd;
xd=x;
nr=n;
int s;
for (s=0;s<15;s++)
{ //
sr=s;
ansd=ansd+pow(-1.0,sr)*factorial(sr+nr)*pow(xd,2.0*sr)/(factorial(sr)*factorial(2.0*s+2.0*n+1));
} //
ansd = ansd * pow(2.0,nr)*pow(xd,nr);
ans=ansd;
return ans;
} //

float THETAx(int n,float x,float y,float z)
{ // x derivative of THETA
float x1,x2,dx=0.01,ans,THETA2=0,THETA1=0;
float r,costheta;
x2=x+0.5*dx; r=sqrt(x2*x2+y*y+z*z);
if (r!=0) { costheta=z/r; THETA2=SphBess1j(n,klguess*r)*LegenP(n,costheta); }
x1=x-0.5*dx; r=sqrt(x1*x1+y*y+z*z);
if (r!=0) { costheta=z/r; THETA1=SphBess1j(n,klguess*r)*LegenP(n,costheta); }
ans=(THETA2-THETA1)/dx;
return ans;
} // x derivative of THETA

float THETAy(int n,float x,float y,float z)
{ // y derivative of THETA
float y1,y2,dx=0.01,ans,THETA2=0,THETA1=0;
float r,costheta;
y2=y+0.5*dx; r=sqrt(x*x+y2*y2+z*z);
if (r!=0) { costheta=z/r; THETA2=SphBess1j(n,klguess*r)*LegenP(n,costheta); }
y1=y-0.5*dx; r=sqrt(x*x+y1*y1+z*z);
if (r!=0) { costheta=z/r; THETA1=SphBess1j(n,klguess*r)*LegenP(n,costheta); }
ans=(THETA2-THETA1)/dx;
return ans;
} // y derivative of THETA

float THETAz(int n,float x,float y,float z)
{ // z derivative of THETA
float z1,z2,dx=0.01,ans,THETA2=0,THETA1=0;
float r,costheta;
z2=z+0.5*dx; r=sqrt(x*x+y*y+z2*z2);
if (r!=0) { costheta=z2/r; THETA2=SphBess1j(n,klguess*r)*LegenP(n,costheta); }
z1=z-0.5*dx; r=sqrt(x*x+y*y+z1*z1);
if (r!=0) { costheta=z1/r; THETA1=SphBess1j(n,klguess*r)*LegenP(n,costheta); }
ans=(THETA2-THETA1)/dx;
return ans;
} // z derivative of THETA

void makescr(void)
{ // makescr
//gotoxy(70,5);cout<<"Make .scr ";
//gotoxy(70,6);cout<<"file? ";
//gotoxy(70,7);cout<<"(y/n/q)  ";
//char ch678;
//ch678=getch();
//if (ch678=='q') {closegraph();exit(0);}
//if (ch678!='y') return;
//char scrfilename[40];
//FILE *scrfout;
//gotoxy(65,8);cout<<"Filename";
//gotoxy(65,9);cout<<"? ";cin>>scrfilename;
int ExtraRowsForTitle=2;
//scrfout=fopen(scrfilename,"wb");
unsigned char ch1b,ch2b,ch3b,ch4b,ch5,ch6;
//int bx0=100;
//int by0=200;
//setcolor(7);
//rectangle(bx0-89-1,by0-89-1,bx0+89+1,by0+89+1);

int scrwid=2*85;
int scrhgt=2*85;
int scrx0=100-85;
int scry0=200-85;
setcolor(15);
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,10);cout<<"saving ";
for (iy=scry0;iy<scry0+scrhgt;iy++)
{ // iy
gotoxy(65,11);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

void mn(char *ch1)
{ // mn
strcpy(matername,ch1);
return;
} // mn

void circleb(float rix,float riy,float r1,float r2,int icol)
{ // circleb
int ix3,iy3,ir3;
float r32,r12,r22;
r12=r1*r1;
r22=r2*r2;
ir3=r2+1.6;
for (ix3=rix-ir3;ix3<=rix+ir3;ix3++)
for (iy3=riy-ir3;iy3<=riy+ir3;iy3++)
{ // ix3 iy3
r32=(ix3-rix)*(ix3-rix)+(iy3-riy)*(iy3-riy);
if ((r32<=r22)&&(r32>=r12))
{ //
putpixel(ix3,iy3,icol);
} //
} // ix3 iy3
return;
} // circleb