//C program to calculate the uncertainty fraction of a scattering system:
#include <stdio.h>
#include <math.h>
#include "random.h"
#include "ndisc.h"
//#define SCATFILE "grscat.init"
#define BMIN 0.
#define BMAX 1.
#define X0 -4.
#define THETA0 0.
#define REPORTFILE "scatrep.txt"
//syntax:
//
//scat_uf mineps maxeps neps minunc maxtrial outfile noisetype [noisepow]
int main(int argc, char *argv[]) {
float b; //impact parameter
float xf, yf; //final x and y values
long nbounce1, nbounce2; //number of bounces
float t; //time delay
float theta; //direction of travel
float epsmin; //minimum epsilon value
float epsmax; //maximum epsilon value
long neps; //number of epsilon values
float eps; //current epsilon value
char *outfile; //output file
FILE *outfs; //output file stream
long minunc; //minimum number of uncertain values
long maxtrial; //maximum number of trials
long nunc; //number of uncertain trials
float f; //uncertainty fraction
char noisetype; //noise type
float std; //noise power
long seed;
Ndisc system; //scattering system
//parse the command line arguments:
if (argc < 9) {
printf("syntax: scatuf mineps maxeps neps minunc maxn outfile noisetype noisepow\n");
printf("\n");
printf("where:\n");
printf(" mineps minimum epsilon value\n");
printf(" maxeps maximum epsilon value\n");
printf(" neps number of values for epsilon (logarithmic progression)\n");
printf(" minunc minimum number of uncertain trials\n");
printf(" maxn maximum number of trials to determine uncertainty fraction\n");
printf(" outfile output file\n");
printf(" noisetype w = white; g = gaussian; n = no noise\n");
printf(" noisepow noise power (standard deviation in radians)\n");
return 1;
}
sscanf(argv[1], "%g", &epsmin);
sscanf(argv[2], "%g", &epsmax);
sscanf(argv[3], "%d", &neps);
sscanf(argv[4], "%d", &minunc);
sscanf(argv[5], "%d", &maxtrial);
outfile=argv[6];
noisetype=argv[7][0];
sscanf(argv[8], "%g", &std);
long j;
//initialise the scattering system:
if (noisetype == 'g') {
system.init(GAUSS);
} else {
system.init(WHITE);
}
//open up the output file:
outfs=fopen(outfile, "w");
if (noisetype == 'n') {
float deps=pow(epsmax/epsmin, 1./(neps-1));
eps=epsmin;
for (long i=0; i<neps; i++) {
nunc=0;
for (j=0; j<maxtrial; j++) {
b=BMIN+ran2(&seed)*(BMAX-BMIN);
xf=X0;
yf=b;
theta=THETA0;
system.scatter(&xf, &yf, &theta, &t, &nbounce1);
b=b+(ran2(&seed)*2-1)*eps;
xf=X0;
yf=b;
theta=THETA0;
system.scatter(&xf, &yf, &theta, &t, &nbounce2);
if (nbounce1 != nbounce2) nunc++;
if (nunc > minunc) break;
}
f=(float) nunc/(float) j;
fprintf(outfs, "%12.8g %12.8g\n", eps, f);
eps*=deps;
}
} else {
float deps=pow(epsmax/epsmin, 1./(neps-1));
eps=epsmin;
for (long i=0; i<neps; i++) {
nunc=0;
for (j=0; j<maxtrial; j++) {
b=BMIN+ran2(&seed)*(BMAX-BMIN);
xf=X0;
yf=b;
theta=THETA0;
system.scatter_noisy(&xf, &yf, &theta, &t, &nbounce1, std);
b=b+(ran2(&seed)*2-1)*eps;
xf=X0;
yf=b;
theta=THETA0;
system.scatter_noisy(&xf, &yf, &theta, &t, &nbounce2, std);
if (nbounce1 != nbounce2) nunc++;
if (nunc > minunc) break;
}
f=(float) nunc/(float) j;
fprintf(outfs, "%12.8f %12.8g\n", eps, f);
eps*=deps;
}
}
fclose(outfs);
}