// 8ptc11.cpp  Aug 2 2002   8 point force
//  This program establishes that:
//           a C11 = ksp1 + 2 ksp2 + 0.66666 k8pt

int Polarization=0; // 0=Long; 1=shear y; 2 = shear z

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

double omega,kx,totime,dt,matom,lambda,CL,uca,Amp,densityrho;
double u(int nx,int ny,int nz,int iaxis);
long int countit;

void main(void)
{ // main
clrscr();
//int gdriver=VGA;
//int gmode=VGAHI;
//initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
//cleardevice();
//setcolor(8);
//rectangle(0,0,639,479);
kx=0.041;
lambda=2.0*3.1415926/kx;
uca=0.59;
densityrho=1.0;
matom=densityrho*uca*uca*uca;
dt=0.1;
omega=2.5;
Amp=0.0001251; // amplitude of wave

int nx,ny,nz;
double sfx,bx0,by0;
sfx=70;
bx0=320;
by0=240;
countit=0;

LBtimestep:
countit++;
totime=totime+dt;
gotoxy(1,1);cout<<"it= "<<countit<<" ";
gotoxy(1,2);cout<<"t= "<<totime<<" s";
nz=0;
for (nx=-3;nx<=3;nx++)
for (ny=-3;ny<=3;ny++)
{ //
double ux,uy,uz;
ux = u(nx,ny,nz,0);
uy = u(nx,ny,nz,1);
uz = u(nx,ny,nz,2);
//putpixel(bx0+sfx*(nx*uca+ux),by0-sfx*(ny*uca+uy),15);
} //
delay(300);
double ksp1,ksp2,k8pt;
double F000x,F000y,F000z,Fmag,r01,r012,dx,dy,dz,Fx,Fy,Fz; // tot
int N2;
ksp1=.0;
ksp2=.30;
k8pt=0.9;
F000x=0;
F000y=0;
F000z=0;
for (nx=-1;nx<=1;nx++)
for (ny=-1;ny<=1;ny++)
for (nz=-1;nz<=1;nz++)
{ //.. nx ny nz
N2=nx*nx+ny*ny+nz*nz;
dx=nx*uca+u(nx,ny,nz,0)-u(0,0,0,0);
dy=ny*uca+u(nx,ny,nz,1)-u(0,0,0,1);
dz=nz*uca+u(nx,ny,nz,2)-u(0,0,0,2);
r012=dx*dx+dy*dy+dz*dz;
r01=sqrt(r012);
Fmag=0;
if (N2==1)
{ // first neighbour
Fmag=ksp1*(r01-uca*sqrt(N2));
} // first neighbour
if (N2==2)
{ // 2nd neigh.
Fmag=ksp2*(r01-uca*sqrt(N2));
} // 2nd
Fx=0;Fy=0;Fz=0;
if (N2!=0)
{ //.
Fx = (dx/r01)*Fmag;
Fy = (dy/r01)*Fmag;
Fz = (dz/r01)*Fmag;
} //.
F000x+=Fx;
F000y+=Fy;
F000z+=Fz;
} //.. nx ny nz
// 8 point force:
//  (on atom 000 from 8 different clusters of 8)
for (nx=-1;nx<=0;nx++)
for (ny=-1;ny<=0;ny++)
for (nz=-1;nz<=0;nz++)
{ //.. nx ny nz
double cm[3]; // Find Center of mass of cluster of 8
cm[0]=0;
cm[1]=0;
cm[2]=0;
int nx1,ny1,nz1;
for (nx1=0;nx1<=1;nx1++)
for (ny1=0;ny1<=1;ny1++)
for (nz1=0;nz1<=1;nz1++)
{ //.. nx1 ny1 nz1
int nx2,ny2,nz2;
nx2=nx+nx1;
ny2=ny+ny1;
nz2=nz+nz1;
cm[0]=cm[0]+uca*nx2 + u(nx2,ny2,nz2,0);
cm[1]=cm[1]+uca*ny2 + u(nx2,ny2,nz2,1);
cm[2]=cm[2]+uca*nz2 + u(nx2,ny2,nz2,2);
} //.. nx1 ny nz1
cm[0]=cm[0]/8.0;
cm[1]=cm[1]/8.0;
cm[2]=cm[2]/8.0;
//putpixel(bx0+sfx*(cm[0]),by0-sfx*(cm[1]),13);
double sum8;
sum8=0;
for (nx1=0;nx1<=1;nx1++)
for (ny1=0;ny1<=1;ny1++)
for (nz1=0;nz1<=1;nz1++)
{ //  nx1 ny1 nz1
int nx2,ny2,nz2;
nx2=nx+nx1;
ny2=ny+ny1;
nz2=nz+nz1;
dx = uca*nx2 + u(nx2,ny2,nz2,0) - cm[0];
dy = uca*ny2 + u(nx2,ny2,nz2,1) - cm[1];
dz = uca*nz2 + u(nx2,ny2,nz2,2) - cm[2];
sum8=sum8+dx*dx+dy*dy+dz*dz;
} //  nx1 ny1 nz1
sum8=sum8-6.0*uca*uca;
gotoxy(30,10);cout<<sum8<<"  ";
sum8=sum8*k8pt/(12.0*uca*uca);
Fx=sum8*(cm[0]-u(0,0,0,0));
Fy=sum8*(cm[1]-u(0,0,0,1));
Fz=sum8*(cm[2]-u(0,0,0,2));
F000x+=Fx;
F000y+=Fy;
F000z+=Fz;
} //.. nx ny nz

gotoxy(60,3);cout<<"Fx "<<F000x<<"  ";
gotoxy(60,4);cout<<"Fy "<<F000y<<"  ";
gotoxy(60,5);cout<<"Fz "<<F000z<<"  ";

double u000;
if (Polarization==0) u000=u(0,0,0,0);
if (Polarization==1) u000=u(0,0,0,1);
if (Polarization==2) u000=u(0,0,0,2);

if (u000!=0)
{ //
double Fxbyx;
if (Polarization==0) Fxbyx=F000x/u(0,0,0,0);
if (Polarization==1) Fxbyx=F000y/u(0,0,0,1);
if (Polarization==2) Fxbyx=F000z/u(0,0,0,2);
gotoxy(60,7);cout<<"Fx/x "<<Fxbyx<<"  ";
double omega2=sqrt(fabs(Fxbyx/matom));
gotoxy(60,8);cout<<"omega "<<omega2<<"  ";
double CL2=omega2/kx;
if (Polarization==0) {gotoxy(60,9);cout<<"CL = "<<CL2<<"  ";}
if (Polarization==1) {gotoxy(60,9);cout<<"CS = "<<CL2<<"  ";}
if (Polarization==2) {gotoxy(60,9);cout<<"CS = "<<CL2<<"   ";}
double C11;
C11 = (CL2*CL2)*densityrho;
if (Polarization==0) {gotoxy(60,10);cout<<"C11 = "<<C11<<"   ";}
if (Polarization==1) {gotoxy(60,10);cout<<"C44 = "<<C11<<"   ";}
if (Polarization==2) {gotoxy(60,10);cout<<"C44 = "<<C11<<"   ";}
if (Polarization==0) {gotoxy(50,11);cout<<"(ksp1+2ksp2)/a = "<<(ksp1+2*ksp2)/uca<<"   ";}
if (Polarization==1) {gotoxy(50,11);cout<<"(ksp2)/a = "<<(ksp2)/uca<<"  ";}
if (Polarization==2) {gotoxy(50,11);cout<<"(ksp2)/a = "<<(ksp2)/uca<<"  ";}
if (k8pt!=0)
{ //
if (Polarization==0) {gotoxy(50,13);cout<<"(??) = "<<(uca*C11-(ksp1+2*ksp2))/k8pt<<"   ";}
} //
// Fx/x = m omega * omega
} //

for (nx=-3;nx<=3;nx++)
for (ny=-3;ny<=3;ny++)
{ //
double ux,uy,uz;
ux = u(nx,ny,nz,0);
uy = u(nx,ny,nz,1);
uz = u(nx,ny,nz,2);
//putpixel(bx0+sfx*(nx*uca+ux),by0-sfx*(ny*uca+uy),0);
} //



if (!kbhit()) goto LBtimestep;
char ch1;
ch1=getch();
if (ch1=='q') return;
//if (ch1=='c') cleardevice();
if (ch1=='k') kx=kx*0.7;
if (ch1=='a') Amp=Amp*0.7;
goto LBtimestep;


getch();
return;
} // main

double u(int nx,int ny,int nz,int iaxis)
{ // u
int LF=Polarization;
double ux;
if ((LF==0)&&(iaxis==0)) return Amp*cos(kx*nx*uca-omega*totime);
if ((LF==0)&&(iaxis==1)) return 0;
if ((LF==0)&&(iaxis==2)) return 0;
if ((LF==1)&&(iaxis==0)) return 0;
if ((LF==1)&&(iaxis==1)) return Amp*cos(kx*nx*uca-omega*totime);
if ((LF==1)&&(iaxis==2)) return 0;
if ((LF==2)&&(iaxis==0)) return 0;
if ((LF==2)&&(iaxis==1)) return 0;
if ((LF==2)&&(iaxis==2)) return Amp*cos(kx*nx*uca-omega*totime);
cout<<" wrong axis/polarization ";
exit(0);
} // u