#include <stdio.h>
#include <stdlib.h>
//#include "e:\files2\sci_lib\graph_to.h"
#include "/home/mills/sci_lib/random.h"
/*#include <math.h>*/
#define CIRCLENUMBER 3
#define PI 3.141592654
/*This include file contains the object definitions for a simulation of the three-disc scattering
system*/
class disc {
float rad;
float cx, cy;
public:
void init(float radius, float centerx, float centery);
void intercept (float begx, float begy, float dirx, float diry,
float *t1, float *t2);
int int_near (float begx, float begy, float dirx, float diry, float *t);
int int_far (float begx, float begy, float dirx, float diry, float *t);
void bounce (float interceptx, float intercepty, float dirx, float diry,
float *reflectxhat, float *reflectyhat);
void plot();
};
/************************************************************************/
class Ndisc {
/*Contains several disc objects which compose the scattering system*/
disc alldiscs[CIRCLENUMBER]; /*pointer to array of disc objects*/
disc boundary; /*defines the outer border at which particle
exits the system*/
float (*noise_func)(long *idum);
long idum; /*seed value for randum number generator*/
int Ndisc::intercept(float x, float y, float dirx, float diry, float *t);
public:
void init(int noise_type);
/*initializes object with 3 discs in the usual pattern*/
/*noise_type: 0=gaussian, 1=white*/
void init2 (int noise_type);
/*initializes object same as experimental apparatus*/
void plot_discs();
/*plots the scattering system*/
#if NOISY
void scatter (float *startx1, float *starty1, float *theta, float *tdelay,
long *bounces, float var);
#else
void scatter (float *startx1, float *starty1, float *theta, float *tdelay,
long *bounces);
#endif
/*simulates a point particle moving through the scattering system*/
};
/************************************************************************/
//graph con; /*Defines graphics conversions*/
/************************************************************************/
void disc::init(float radius, float centerx, float centery) {
rad=radius;
cx=centerx;
cy=centery;
}
/*______________________________________________________________________*/
void disc::intercept(float begx, float begy, float dirx, float diry,
float *t1, float *t2)
/*Given the origin of a line, its direction vector, the center of a circle,
its radius vector, find the intercept points between that line and the
circle, and return in the form of the parameters of the line, t1 & t2.
Returns -1 for each parameter if there are no intercepts*/
{
float X; /*intermediate value in calculation*/
float Y; /* " */
float B; /* " */
float Atimes2; /* " */
float discriminant; /*Square root of the discriminant in quadratic*/
X=begx-cx;
Y=begy-cy;
B=dirx*X+diry*Y;
Atimes2=2*(dirx*dirx+diry*diry);
discriminant=-pow(dirx*Y-diry*X,2)+pow(dirx*rad,2)+pow(diry*rad,2);
if (discriminant >= 0) {
*t1=(-2*B+2*sqrt(discriminant))/Atimes2;
*t2=(-2*B-2*sqrt(discriminant))/Atimes2;
} else {
*t1=-1;
*t2=-1;
}
}
/*______________________________________________________________________*/
int disc::int_near(float begx, float begy, float dirx, float diry, float *t)
/*Given the origin of a line, its direction vector, the center of a circle,
its radius vector, returns the closest intercept in the form of the
parameter of the line, t.
Returns -1 for each parameter if there are no intercepts*/
{
float X; /*intermediate value in calculation*/
float Y; /* " */
float B; /* " */
float Atimes2; /* " */
float discriminant; /*Square root of the discriminant in quadratic*/
X=begx-cx;
Y=begy-cy;
B=dirx*X+diry*Y;
Atimes2=2*(dirx*dirx+diry*diry);
discriminant=-pow(dirx*Y-diry*X,2)+pow(dirx*rad,2)+pow(diry*rad,2);
if (discriminant >= 0) {
*t=(-2*B-2*sqrt(discriminant))/Atimes2;
return 1;
} else {
*t=-1;
return 0;
}
}
/*______________________________________________________________________*/
int disc::int_far(float begx, float begy, float dirx, float diry, float *t)
/*Given the origin of a line, its direction vector, the center of a circle,
its radius vector, returns the farthest intercept in the form of the
parameter of the line, t.
/*Returns 0, *t=-1 if there are no intercepts 1 if there are*/
{
float X; /*intermediate value in calculation*/
float Y; /* " */
float B; /* " */
float Atimes2; /* " */
float discriminant; /*Square root of the discriminant in quadratic*/
X=begx-cx;
Y=begy-cy;
B=dirx*X+diry*Y;
Atimes2=2*(dirx*dirx+diry*diry);
discriminant=-pow(dirx*Y-diry*X,2)+pow(dirx*rad,2)+pow(diry*rad,2);
if (discriminant >= 0) {
*t=(-2*B+2*sqrt(discriminant))/Atimes2;
return 1;
} else {
*t=-1;
return 0;
}
}
/*______________________________________________________________________*/
void disc::bounce(float interceptx, float intercepty, float dirx, float diry,
float *reflectxhat, float *reflectyhat) {
float normx, normy; /*x,y coords of normal to radius vector*/
float reflectlength; /*projection of direction vector on normal*/
normx=intercepty-cy;
normy=cx-interceptx;
reflectlength=(normx*dirx+normy*diry)/(normx*normx+normy*normy);
*reflectxhat=-dirx+2*reflectlength*normx;
*reflectyhat=-diry+2*reflectlength*normy;
}
void disc::plot() {
// con.disc(cx,cy,rad);
}
/**************************************************************************/
void Ndisc::init (int noise_type) {
if (noise_type==1) {
noise_func=&white;
} else {
noise_func=&gasdev;
}
idum=-784521258;
alldiscs[0].init(1,-0.7216878,1.25);
alldiscs[1].init(1,1.4433756,0);
alldiscs[2].init(1,-0.7216878,-1.25);
boundary.init(5,0,0);
}
/*______________________________________________________________________*/
void Ndisc::init2 (int noise_type) {
if (noise_type==1) {
noise_func=&white;
} else {
noise_func=&gasdev;
}
idum=-784521258;
alldiscs[0].init(2.5,-1.847520861,3.2);
alldiscs[1].init(2.5,3.69041723,0);
alldiscs[2].init(2.5,-1.847520861,-3.2);
boundary.init(7,0,0);
}
/*______________________________________________________________________*/
void Ndisc::plot_discs() {
int i;
/* con.window(-3,-2.5,3,2.5);*/
// con.window(-4,-3,4,3);
/* con.window(-9,-7,9,7);*/
/* con.pos(-3,0);
// con.linerel(3,0);
// con.pos(0,-2.5);
// con.linerel(0,2.5);*/
for (i=0;i<=CIRCLENUMBER-1;i++) {
alldiscs[i].plot();
}
}
/*________________________________________________________________________*/
int Ndisc::intercept(float x, float y, float dirx, float diry, float *t) {
/*Calculates the distance to the nearest point where the particle intercepts
a disc. Returns this value in *t. Returns the number of the circle. If
this is greater than CIRCLENUMBER then the particle will exit the system.*/
float ttest; /*test value for distance to intercept*/
float tdum; /*value for intercept which is thrown away*/
int whichcircle; /*the disc which the particle intercepts*/
int a; /*counts through the discs*/
int exit_system; /*particle exits the system*/
exit_system=1;
for (a=0;a<=CIRCLENUMBER-1;a++) {
/* alldiscs[a].intercept (x,y,dirx,diry,&tdum,&ttest);*/
if (alldiscs[a].int_near (x,y,dirx,diry,&ttest)&&(ttest>0.0)) {
if (exit_system) {
*t=ttest;
exit_system=0;
whichcircle=a;
} else if ((*t>ttest)&&(ttest>0.0)) {
*t=ttest;
whichcircle=a;
}
}
}
if (exit_system) {
/* boundary.intercept(x,y,dirx,diry,t,&tdum);*/
if (!boundary.int_far (x,y,dirx,diry,t)) {
fprintf(stderr, "Improper initial conditions");
exit(1);
}
whichcircle=CIRCLENUMBER;
}
return whichcircle;
}
/*______________________________________________________________________*/
void Ndisc::scatter (float *startx, float *starty, float *theta, float *tdelay,
long *bounces
#if NOISY
, float stddev)
#else
)
#endif
/*Given starting coords and angle, find time for particle to exit system
final exit pint is returned in startx1, starty1 while angle at which the
particle exits the system is returned in theta.*/
{
float tmin; /*final value for distance to intercept*/
float interceptx; /*coordinate of interception with disc*/
float intercepty;
float dirxhat; /*momentum (normalized)*/
float diryhat;
float deltatheta; /*random change in reflection angle for noise*/
float cosdel; /*cos of deltatheta*/
float sindel; /*sin " " */
float theta1; /*exit angle*/
int whichcircle;
int a;
float reflectxhat, reflectyhat; /*momentum after reflection*/
dirxhat=cos(*theta);
diryhat=sin(*theta);
*bounces=0;
*tdelay=0;
whichcircle=0;
// con.pos(*startx,*starty);
while (whichcircle!=CIRCLENUMBER) {
whichcircle=intercept(*startx,*starty,dirxhat,diryhat,&tmin);
if (whichcircle==CIRCLENUMBER) {
*startx=(*startx+dirxhat*tmin);
*starty=(*starty+diryhat*tmin);
(*tdelay)+=tmin;
// con.linerel(*startx, *starty);
} else {
(*tdelay)+=tmin;
(*bounces)++;
/* printf ("%5i\n",*bounces);*/
interceptx=(*startx+dirxhat*tmin);
intercepty=(*starty+diryhat*tmin);
alldiscs[whichcircle].bounce(interceptx,intercepty,dirxhat,diryhat,
&reflectxhat,&reflectyhat);
#if NOISY
deltatheta=stddev*(*noise_func)(&idum);
sindel=sin(deltatheta);
cosdel=cos(deltatheta);
dirxhat=reflectxhat*cosdel-reflectyhat*sindel;
diryhat=reflectyhat*sindel+reflectyhat*cosdel;
#else
dirxhat=reflectxhat;
diryhat=reflectyhat;
#endif
*startx=interceptx;
*starty=intercepty;
// con.linerel(*startx,*starty);
}
}
theta1=atan(diryhat/dirxhat);
if (dirxhat<0.0) {*theta=theta1-PI*fabs(theta1)/theta1;}
else {*theta=theta1;}
}