-
Notifications
You must be signed in to change notification settings - Fork 2
/
Source.cpp
88 lines (71 loc) · 1.96 KB
/
Source.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include <cmath>
#include <vector>
#include <tuple>
#include "Random.h"
#include "Source.h"
#include "Particle.h"
unsigned int Source::groupSample(std::vector<double> groupProbability)
{
if(groupProbability.size() > 1) {
double rand = Urand();
double c = 0;
for(unsigned int i = 0; i < groupProbability.size(); i++)
{
c += groupProbability[i];
if(c > rand)
{
return (i+1);
}
}
}
return(1);
}
Part_ptr setSourcePoint::sample(){
double pi = acos(-1.);
auto group = groupSample(groupProbability);
point pos = point(x0,y0,z0);
double mu = 2 * Urand() - 1;
double phi = 2 * pi*Urand();
double omegaX=mu;
double omegaY=sin(acos(mu))*cos(phi);
double omegaZ=sin(acos(mu))*sin(phi);
point dir = point(omegaX,omegaY,omegaZ);
Part_ptr p = std::make_shared<Particle>(pos, dir, group );
return p;
}
Part_ptr setSourceSphere::sample(){
double pi = acos(-1.);
//Radius of the new particle
//double radius = pow((pow(radInner,3.0) + Urand()*(pow(radOuter,3.0)-pow(radInner,3.0))),(1. / 3.));
//double mu = 2.0 * Urand() - 1.0;
//double phi = 2.0 * pi*Urand();
//double x=radius*sqrt(1-pow(mu,2.))*cos(phi)+x0;
//double y=radius*sqrt(1-pow(mu,2.))*sin(phi)+y0;
//double z=radius*mu+z0;
//std::cout << "mu: " << mu << " x: " << x << " y: " << y << " z: " << z << std::endl;
//std::cout << std::endl;
double x;
double y;
double z;
bool reject = true;
while(reject)
{
x = 2*Urand()*radOuter;
y = 2*Urand()*radOuter;
z = 2*Urand()*radOuter;
double dist = sqrt(x*x+y*y+z*z);
if(dist < radOuter)
reject = false;
}
auto group = groupSample(groupProbability);
point pos = point(x,y,z);
// direction sampling
double mu = 2 * Urand() - 1;
double phi = 2 * pi*Urand();
double omegaX=mu;
double omegaY=sin(acos(mu))*cos(phi);
double omegaZ=sin(acos(mu))*sin(phi);
point dir = point(omegaX,omegaY,omegaZ);
Part_ptr p = std::make_shared<Particle>(pos, dir, group );
return p;
}