#define VERN "herve15c.cpp"
//  15c
// herve15b   add d electron susceptibility term
// herve15.cpp  May 9, 2003
//  - try to reproduce figure 1.5 in Herve portales' thesis.  these.pdf

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


void main(void)
{ // main
int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();
setcolor(10);
rectangle(0,0,639,479);
float hervexi;
float semiaxisa,semiaxisb,semiaxisc;
float ecc;

int ipx,ipy1,ipy2,ipy3,ipy4;
int bx0=150,by0=300;
float sfx=40.0,sfy=50,ymin=0.5,ymax=4.0;

setcolor(10);
line(bx0+sfx*log(0.088),by0-sfy*ymin,bx0+sfx*log(10.0),by0-sfy*ymin);
line(bx0+sfx*log(0.088),by0-sfy*ymin,bx0+sfx*log(0.088),by0-sfy*ymax);

for (hervexi=0.125;hervexi<=8.0;hervexi=hervexi*2.0)
{ // hervexi
ipx = bx0 + sfx*log(hervexi);
setcolor(10);
line(ipx,by0-sfy*ymin,ipx,by0+5-sfy*ymin);
setcolor(8);
line(ipx,by0-sfy*ymin,ipx,by0+-sfy*ymax);
} // hervexi

setcolor(15);
settextstyle(2,0,5);
outtextxy(-2+bx0+sfx*log(1.0),by0-17,"1");
outtextxy(-2+bx0+sfx*log(2.0),by0-17,"2");
outtextxy(-2+bx0+sfx*log(4.0),by0-17,"4");
outtextxy(-2+bx0+sfx*log(8.0),by0-17,"8");
outtextxy(-5+bx0+sfx*log(0.5),by0-17,".5");
outtextxy(-8+bx0+sfx*log(0.25),by0-17,".25");
outtextxy(-11+bx0+sfx*log(0.125),by0-17,".125");

hervexi=1.0;
ipx = bx0 + sfx*log(hervexi);
setcolor(9);
line(ipx,by0-sfy*ymin,ipx,by0-sfy*ymax);
hervexi=10.0;
ipx = bx0 + sfx*log(hervexi);
setcolor(10);
line(ipx,by0-sfy*ymin,ipx,by0-sfy*ymax);


float omega1;
gotoxy(2,3);cout<<VERN<<"   Ag ellipsoid plasmon";
gotoxy(2,4);cout<<"After Figure 1.3 on page 24 of thesis";
gotoxy(2,5);cout<<"by Herv‚ Portales";


setcolor(15);
settextstyle(2,0,5);
outtextxy(-15+bx0+sfx*log(0.088),-8+by0-sfy*1.0,"1");
outtextxy(-15+bx0+sfx*log(0.088),-8+by0-sfy*2.0,"2");
outtextxy(-15+bx0+sfx*log(0.088),-8+by0-sfy*3.0,"3");
outtextxy(-15+bx0+sfx*log(0.088),-8+by0-sfy*4.0,"4");

for (omega1=0.5;omega1<=4.0;omega1=omega1+0.5)
{ // omega1
setcolor(10);
line(bx0+sfx*log(0.088),by0-sfy*omega1,-4+bx0+sfx*log(0.088),by0-sfy*omega1);
} // omega1

for (omega1=1.0;omega1<=4.0;omega1=omega1+1.0)
{ // omega1
setcolor(10);
line(bx0+sfx*log(0.088),by0-sfy*omega1,-7+bx0+sfx*log(0.088),by0-sfy*omega1);
setcolor(8);
line(bx0+sfx*log(0.088),by0-sfy*omega1,bx0+sfx*log(10),by0-sfy*omega1);
} // omega1

gotoxy(7,20);cout<<"prolate";
gotoxy(25,20);cout<<"oblate";
gotoxy(2,11);cout<<"freq";
gotoxy(3,12);cout<<"eV";

for (hervexi=0.088;hervexi<=10.0;hervexi=hervexi*1.011)
{ // hervexi
semiaxisb=1;
semiaxisc=1;
semiaxisa=semiaxisb/hervexi;
gotoxy(40,1);cout<<"xi = "<<hervexi<<"  ";
gotoxy(40,2);cout<<"a = "<<semiaxisa<<"  ";
gotoxy(40,3);cout<<"b = "<<semiaxisb<<"  ";
gotoxy(40,4);cout<<"c = "<<semiaxisc<<"  ";
float dpfnx,dpfny,dpfnz;

if (semiaxisa>=semiaxisb)
{ // prolate
gotoxy(40,5);cout<<"Prolate spheroid";
ecc = sqrt(1.0-(semiaxisb/semiaxisa)*(semiaxisb/semiaxisa));
gotoxy(40,6);cout<<"ecc = "<<ecc<<" ";
float e1=ecc,e2=ecc*ecc,e3=ecc*ecc*ecc;
dpfnx = ((1-e2)/(2*e3))*(log((1+e1)/(1-e1))-2*e1);
gotoxy(40,7);cout<<"n(x) = "<<dpfnx<<" ";
dpfny = 0.5*(1-dpfnx);
dpfnz = 0.5*(1-dpfnx);
gotoxy(40,8);cout<<"n(y) = "<<dpfny<<" ";
gotoxy(40,9);cout<<"n(z) = "<<dpfnz<<" ";
} // prolate
if (semiaxisa<semiaxisb)
{ // oblate
gotoxy(40,5);cout<<"Oblate spheroid";
ecc = sqrt((semiaxisc/semiaxisa)*(semiaxisc/semiaxisa)-1.0);
gotoxy(40,6);cout<<"ecc = "<<ecc<<" ";
float e1=ecc,e2=ecc*ecc,e3=ecc*ecc*ecc;
dpfnx = ((1+e2)/e3)*(e1-atan(e1));
gotoxy(40,7);cout<<"n(x) = "<<dpfnx<<" ";
dpfny = 0.5*(1-dpfnx);
dpfnz = 0.5*(1-dpfnx);
gotoxy(40,8);cout<<"n(y) = "<<dpfny<<" ";
gotoxy(40,9);cout<<"n(z) = "<<dpfnz<<" ";
} // oblate

float omegap_eV = 8.98; // plasma frequency in electron volts
float omega_eV_para;
float omega_eV_perp;
float omega_eV_para_al2o3;
float omega_eV_perp_al2o3;
float epsrm = 1.0;  // relative permittivity of matrix material
float epsrm_al2o3 = 3.0;  // relative permittivity of Al2O3 material
float xid = 3.0; // d electron susceptibility
// See figure 1.3 in Herve Portales' thesis to see actual frequency
// dependence of d electron contribution to susceptibility.
omega_eV_para = omegap_eV / sqrt(xid+1.0-epsrm*(1.0-1.0/dpfnx));
omega_eV_perp = omegap_eV / sqrt(xid+1.0-epsrm*(1.0-1.0/dpfnz));
omega_eV_para_al2o3 = omegap_eV / sqrt(xid+1.0-epsrm_al2o3*(1.0-1.0/dpfnx));
omega_eV_perp_al2o3 = omegap_eV / sqrt(xid+1.0-epsrm_al2o3*(1.0-1.0/dpfnz));
//gotoxy(40,11);cout<<"para: omega = "<<omega_eV_para<<" eV ";
//gotoxy(40,12);cout<<"perp: omega = "<<omega_eV_perp<<" eV ";
ipx = bx0 + sfx*log(hervexi);
ipy1 = by0 - sfy*omega_eV_para;
ipy2 = by0 - sfy*omega_eV_perp;
ipy3 = by0 - sfy*omega_eV_para_al2o3;
ipy4 = by0 - sfy*omega_eV_perp_al2o3;
putpixel(ipx,ipy1,14);
putpixel(ipx,ipy2,11);
putpixel(ipx,ipy3,13);
putpixel(ipx,ipy4,12);
} // hervexi

setcolor(15);
outtextxy(250,100-33,"vacuum");
outtextxy(250,100-18,"parall");
setcolor(14);
line(250,100,290,100);

setcolor(15);
outtextxy(250,140-33,"Al2O3 ");
outtextxy(250,140-18,"parall");
setcolor(13);
line(250,140,290,140);

setcolor(15);
outtextxy(250,230-33,"vacuum");
outtextxy(250,230-18,"perpen");
setcolor(11);
line(250,230,290,230);

setcolor(15);
outtextxy(250,270-33,"Al2O3");
outtextxy(250,270-18,"perpen");
setcolor(12);
line(250,270,290,270);

makescr();
gotoxy(52,25);cout<<"done....";
//getch();
return;
} // main