// finds.cpp  Aug 6 2002
//  Find complex vibration frequency of an isotropic elastic sphere
//  surrounded by a fluid.

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

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

int gdriver=VGA;
int gmode=VGAHI;
initgraph(&gdriver,&gmode,"c:\\tc\\bgi");
cleardevice();
complex cs,cf1,cf2,czai,ci,ckout,ckr;
float acs;
int ix,iy,icol;
float sfx=30;
float CLa,CTa,rhoa,rhob,Cout;
CLa=3830;
CTa=2640;
rhoa=5680;  // InAs
float Rsph=0.5*55*1e-10; // Radius of sphere
//rhob=860;  // toluene
//Cout=1100;  // toluene (??)
rhob=1110;  // polyvinyl butyral
Cout=2350;  //   polyvinyl butyral

setcolor(15);
float rr,ri,rrbest,ribest,bestabs;
bestabs=1e10;
for (rr=2.1;rr<2.3;rr=rr+0.001)
for (ri=-.3;ri<0;ri=ri+0.001)
{ //
cs = complex(rr,ri);
cf1=tan(cs)/cs;
ci=complex(0,1.0);
ckout=cs*CLa/(Cout*Rsph);
ckr=ckout*Rsph;
czai=rhob*Cout*(ckr*ckr-ci*ckr)/(1.0+(ckr*ckr));
//czai=0;
cf2=(4.0*(CTa*CTa/(CLa*CLa))+ci*(czai*cs/(rhoa*CLa)))*(tan(cs)-cs)/(cs*cs*cs);
acs = 15*abs(cf1-cf2);
if (acs<bestabs) {
bestabs=acs;
rrbest=rr;
ribest=ri;
}
icol=acs;
if (icol>15) icol=15;
setfillstyle(SOLID_FILL,icol);
int ix1=220+sfx*rr;
int iy1=140-sfx*ri;
bar(ix1-1,iy1-1,ix1+1,iy1+1);
} //

cout.precision(3);
gotoxy(5,18);cout<<"Best s = "<<rrbest<<" "<<ribest<<"  ";
float fHz=(CLa*rrbest/Rsph*2.0*3.14159);
gotoxy(5,20);cout<<"Real f = "<<fHz<<" Hz ";
gotoxy(5,21);cout<<"Damping time = "<<-Rsph/(CLa*ribest)<<" s ";

getch();
} // main