7 static const double G = 6.67384e-11;
8 static const double SM = 1.9891e30;
9 static const double PC = 3.08567758e16;
10 static const double C = 299792458;
11 static const double Pi = 3.14159265358979312;
28 double normm = (beta-1)/(1./pow(
minM, beta-1) - 1./pow(
maxM, beta-1));
38 double dr = (rmax - rmin)/nBins;
39 double r = rmin + dr/2;
42 double const bmaxc = pow(340*
Pi/3, 1./7);
44 dGammadmdMdr =
new TH3F(
"d",
"", nBins, mmin, mmax, nBins, mmin, mmax, nBins, rmin, rmax);
48 double vofr = sqrt(mSMBH/rr);
51 for(
int j=0;
j<nBins-1; ++
j){
52 double m1 = mmin + (
j+0.5)*dm;
54 double norm1 = tmp/( pow(rrmax, tmp) - pow(rrmin, tmp) )/4/
Pi ;
57 double m2 = mmin + (
k+0.5)*dm;
59 double eta = m1*m2/M/
M;
61 double bmax = bmaxc*M*pow(eta, 1./7)/pow(vofr, 9./7);
62 double sigmaC =
Pi*bmax*bmax;
64 tmp = 3./2 - p0*m2/
maxM;
65 double norm2 = tmp/( pow(rrmax, tmp) - pow(rrmin, tmp) )/4/
Pi ;
66 double bhDensityProd = normm*normm*norm1*norm2/pow(rr, 3 + p0*M/
maxM)/pow(m1*m2, beta);
68 double rate = bhDensityProd*sigmaC*vofr*4*
Pi*rr*rr;
98 static double g(
double e)
99 {
const double c1 = 12./19;
100 const double c2 = 121./304;
101 const double c3 = 870./2299;
103 return exp(c1*
log(e)+c3*
log(1+c2*e2))/(1-e2);
110 static TF1
F(
"F",
"exp(1.52631578947368429*log(x) + 0.513701609395389336*log(1+0.398026315789473673*x*x) - 1.5*log(1-x*x))", 0, 1);
112 double tM = (m1+
m2)*
SM;
113 double rM = m1*m2/(m1+
m2)*
SM;
114 double e = (ra-
rp)/(ra+rp);
118 double aux =
C*a/
g(e)/
G/sqrt(tM);
121 return F.Integral(0, e)*aux*15./304*
C*
G/rM;
126 {
double tM = 2.8, rM = 0.7;
127 double T0 = 7.75*3600;
129 double a0 = exp(
log(
G*tM*
SM/(4*
Pi*
Pi/T0/T0) )/3 );
132 double rp = (1-
e)*a0;
133 double ra = 2*a0 -
rp;
139 double dade = ((12*aa*(1. + (73./24)*ee*ee + (37./96)*pow(ee,4.)))/(19.*ee*(1-ee*ee)*(1. + (121./304)*ee*ee)));
146 double a = (ra +
rp)/2.;
147 double e = (ra-
rp)/(ra + rp);
148 double totMass = m1+
m2;
150 double k1, k2, k3, k4, k5,
K=1.e10;
152 while(
fabs(1.-k1/K)>.01 ||
fabs(1.-k2/K)>.01 ||
fabs(1.-k3/K)>.01 ||
fabs(1-k4/K)>.01 ||
fabs(1.-k5/K)>.01) {
154 k2 = de *
DiffEq(e + de/3., a + k1/3.);
155 k3 = de *
DiffEq(e + de/3., a + k1/6. + k2/6.);
156 k4 = de *
DiffEq(e + de/2., a + k1/8. + 3.*k3/8.);
157 k5 = de *
DiffEq(e + de, a - k1/2. -3*k3/2. + 2*k4);
158 da = (k1 + 4.*k4 + k5)/6.;
160 K = (k1 + k2 + k3 +k4 +k5)/5.;
165 k2 = de *
DiffEq(e + de/3., a + k1/3.);
166 k3 = de *
DiffEq(e + de/3., a + k1/6. + k2/6.);
167 k4 = de *
DiffEq(e + de/2., a + k1/8. + 3.*k3/8.);
168 k5 = de *
DiffEq(e + de, a - k1/2. -3*k3/2. + 2*k4);
169 da = (k1 + 4.*k4 + k5)/6.;
174 if(da/a>.01) de /=2.;
175 if(da/a<.001) de *=2.;
177 double freq = sqrt(0.5/rp)/(rp*
Pi*totMass*
TimeUnit);
184 {
const double dE_const = 85*
Pi/12/sqrt(2);
186 double M ,
m, eta, rr, vofr, E0, rpm;
196 vofr = sqrt(
smbhM/rr);
199 rpm = pow(dE_const*eta*m/E0, 2./7)*
M;
203 double rndnr =
rnd.Uniform();
209 double L0 = m*sqrt(2*rp*M + pow(vofr*rp, 2));
211 double deltaE = -dE_const*eta*m*pow(M/rp, 3.5);
212 double deltaL = -6*
Pi*pow(M*m/rp, 2);
215 double lom = (L0 + deltaL)/m;
217 double aux2 = sqrt( aux*aux + 2*(E0+deltaE)/m );
219 ra = lom/(aux - aux2);
221 rp = lom/(aux + aux2);
237 if(fn)f=fopen(fn,
"w");
238 for(
int i=0;
i<
n; ++
i){
244 fprintf(f,
"%d %lf %lf %lf %.8lf\n",
i, m1, m2, rp, e);
static double g(double e)
fprintf(stdout,"start=%f duration=%f rate=%f\n", x.start(), x.size()/x.rate(), x.rate())
wavearray< double > a(hp.size())
cout<< endl;cout<< "ts size = "<< ts.size()<< " ts rate = "<< ts.rate()<< endl;tf.Forward(ts, wdm);int levels=tf.getLevel();cout<< "tf size = "<< tf.size()<< endl;double dF=tf.resolution();double dT=1./(2 *dF);cout<< "rate(hz) : "<< RATE<< "\t layers : "<< nLAYERS<< "\t dF(hz) : "<< dF<< "\t dT(ms) : "<< dT *1000.<< endl;int itime=TIME_PIXEL_INDEX;int ifreq=FREQ_PIXEL_INDEX;int index=(levels+1)*itime+ifreq;double time=itime *dT;double freq=(ifreq >0)?ifreq *dF:dF/4;cout<< endl;cout<< "PIXEL TIME = "<< time<< " sec "<< endl;cout<< "PIXEL FREQ = "<< freq<< " Hz "<< endl;cout<< endl;wavearray< double > x
void generateEvents(int n, char *fn=0)
static const double DistUnit
void setFreqCutoff(double f)
static const double TimeUnit
static double Duration(double m1, double m2, double rp, double ra)
GNGen(double mSMBH, double mmin, double mmax, double beta=2)
double generateEvent(double &m1, double &m2, double &rp, double &e)
double fabs(const Complex &x)
double DiffEq(double ee, double aa)
void EvolveRa(double m1, double m2, double &rp, double &ra)