8 static double r,
pr,
phi,
tau,
r2,
a2,
r3,
Delta,
R02,
Q,
A,
B,
r_dot,
w,
pr_dot,
Edot,
Ldot;
16 double& r_dotdot,
double& r_dotdotdot,
double& wdot,
double& wdotdot)
27 double Gdot = E*(3*r2+
a2)*
r_dot;
28 double Qdot = (Fdot*G-F*Gdot)/G/G;
30 double Hdot = L*(
Q*
r_dot+Qdot*
r);
33 wdot = (Hdot*I-H*Idot)/I/I;
34 double Adot = 2*
w*wdot*(
r3-
a2)+
w*
w*(3*r2*
r_dot)+2*a*wdot;
35 double Bdot = (-3*
a2/r2/r2+2/
r3)*r_dot;
38 double Cdot = 4*(r-1)*(r2+
a2-2*r)*
r_dot;
39 double D = E*
R02*r2-2*a*L*
r;
40 double Ddot = E*(4*
r3+2*
a2*
r)*r_dot+r_dot*(2*
a2*E-2*a*L);
42 double alphadot = (Cdot*D-C*Ddot)/D/D;
43 r_dotdot = alphadot*
pr + alpha*
pr_dot;
44 double Cdotdot = 4*(r-1)*(r2+
a2-2*r)*r_dotdot + 4*(
r_dot)*(r2+
a2-2*r)*r_dot
45 + 8*(r-1)*(r-1)*r_dot*
r_dot;
46 double Ddotdot = E*(4*
r3+2*
a2*
r)*r_dotdot+ E*(12*r2+2*
a2)*r_dot*r_dot +
47 r_dotdot*(2*
a2*E-2*a*
L);
48 double alphadotdot = (Cdotdot*D-C*Ddotdot)/D/D - 2*(Cdot*D-C*Ddot)*Ddot/D/D/
D;
49 double pr_dotdot = Adot/r2/
Q-
A*(2*r*r_dot*
Q+r2*Qdot)/(r2*
Q)/(r2*
Q)+Bdot*
pr*
pr*Q
51 r_dotdotdot = alphadotdot*
pr + 2*alphadot*pr_dot + alpha*pr_dotdot;
52 double Fdotdot = (3*r2+
a2-4*
r)*r_dotdot+(6*r-4)*r_dot*
r_dot;
53 double Gdotdot = E*(3*r2+
a2)*r_dotdot + E*6*r*r_dot*r_dot;
54 double Qdotdot = (Fdotdot*G-F*Gdotdot)/G/G-2*(Fdot*G-F*Gdot)*Gdot/G/G/
G;
55 double Hdotdot = L*(Q*r_dotdot+Qdotdot*r+2*Qdot*
r_dot);
56 double Idotdot = (3*r2+
a2)*r_dotdot + 6*r*r_dot*r_dot;
57 wdotdot = (Hdotdot*I-H*Idotdot)/I/I-2*(Hdot*I-H*Idot)*Idot/I/I/I;
62 double phi_t,
double phi_tt,
double phi_ttt,
double (&I_tt)[3][3],
double (&I_ttt)[3][3])
65 double X_tt = 2*r_t*r_t+2*r*r_tt;
66 double X_ttt = 6*r_t*r_tt+2*r*r_ttt;
68 double cos2phi = cos(2*phi);
69 double sin2phi = sin(2*phi);
72 double A_t = -2*sin2phi*phi_t;
73 double A_tt = -4*cos2phi*phi_t*phi_t-2*sin2phi*phi_tt;
74 double A_ttt = 8*sin2phi*pow(phi_t,3)-8*cos2phi*phi_t*phi_tt
75 -4*cos2phi*phi_t*phi_tt-2*sin2phi*phi_ttt;
78 double B_t = 2*cos2phi*phi_t;
79 double B_tt = -4*sin2phi*phi_t*phi_t+2*cos2phi*phi_tt;
80 double B_ttt = -8*cos2phi*pow(phi_t,3) - 8*sin2phi*phi_t*phi_tt
81 -4*sin2phi*phi_t*phi_tt + 2*cos2phi*phi_ttt;
83 double M[3][3] = {1+3*
A,3*
B,0, 3*
B,1-3*
A,0, 0,0,-2};
84 double M_t[3][3] = {3*A_t, 3*B_t, 0, 3*B_t, -3*A_t, 0, 0,0,0};
85 double M_tt[3][3] = {3*A_tt, 3*B_tt, 0, 3*B_tt, -3*A_tt, 0, 0,0,0};
86 double M_ttt[3][3] ={3*A_ttt, 3*B_ttt, 0, 3*B_ttt, -3*A_ttt, 0, 0,0,0};
88 for(
int i=0;
i<3; ++
i)
for(
int j=0;
j<3; ++
j){
89 I_tt[
i][
j] = (X_tt*M[
i][
j]+X*M_tt[
i][
j]+2*X_t*M_t[
i][
j])/6;
90 I_ttt[
i][
j] = (X_ttt*M[
i][
j]+X*M_ttt[
i][
j]+3*X_tt*M_t[
i][
j]+3*X_t*M_tt[
i][
j])/6;
94 void quad_loss(
double r,
double r_t,
double r_tt,
double r_ttt,
double phi,
95 double phi_t,
double phi_tt,
double phi_ttt,
double& E_dot,
double& Lz_dot)
105 {
double I_tt[3][3], I_ttt[3][3];
109 for(
int i=0;
i<3; ++
i)
for(
int j=0;
j<3; ++
j) sum +=I_ttt[
i][
j]*I_ttt[
i][
j];
115 for(
int i=0;
i<3; ++
i)sum += I_tt[0][
i]*I_ttt[
i][1]- I_tt[1][
i]*I_ttt[
i][0];
120 void f_util(
double*
y,
double Jspin,
double mu)
127 double a = Jspin +mu*
L;
144 double r_dotdot, r_dotdotdot, wdot, wdotdot;
155 void h_quad(
double r,
double r_t,
double r_tt,
double r_ttt,
156 double phi,
double phi_t,
double phi_tt,
double phi_ttt,
157 double& hplus,
double& hcross)
161 reduced_quad_dot(r, r_t, r_tt, r_ttt, phi, phi_t, phi_tt, phi_ttt, I_tt, I_ttt);
163 hplus = 2*I_tt[0][0];
164 hcross = 2*I_tt[0][1];
169 int f(realtype
t, N_Vector
y, N_Vector ydot,
void* ud)
170 { realtype *ydata, *ydotdata;
171 realtype* udr = (realtype*)ud;
172 ydata = NV_DATA_S(y);
173 ydotdata = NV_DATA_S(ydot);
175 f_util(ydata, udr[0], udr[1]);
186 double mu,
double dt)
187 {
if(dir!=1 && dir!=-1){
188 printf(
"Error... dir must be +-1\n");
191 y = N_VNew_Serial(6);
206 cvode_mem =CVodeCreate(CV_ADAMS, CV_FUNCTIONAL);
208 printf(
"cvode_mem initialization failed\n");
214 printf(
"cvodeInit failed\n");
218 if(CVodeSStolerances(
cvode_mem, 1
e-4, 1
e-7)!= CV_SUCCESS){
219 printf(
"failed setting tolerances\n");
240 {
return x<y ? x :
y;}
244 this->r =
new double[N+1];
245 this->
pr =
new double[N+1];
246 this->phi =
new double[N+1];
247 this->
tau =
new double[N+1];
248 this->
Et =
new double[N+1];
249 this->
Lt =
new double[N+1];
250 this->
t =
new double[N+1];
252 this->r[0] =
ydata[0];
254 this->phi[0] =
ydata[2];
261 if(stop_alr) r_stop = 2*( 1+cos( 2*(acos(-
min(
udata[0] +
udata[1]*this->
Lt[0], 1))/3 )));
262 if(rmin>r_stop)r_stop = rmin;
264 while(i<N && r[i]>r_stop){
275 r_stop = 2*( 1+cos( 2*(acos(-
min(
udata[0] +
udata[1]*this->Lt[i], 1))/3 )));
276 if(rmin>r_stop)r_stop = rmin;
284 double r_tt, r_ttt, phi_tt, phi_ttt;
287 for(i = 0; i<
N; ++
i){
305 h_quad(::r,
r_dot, r_tt, r_ttt, ::phi,
w, phi_tt, phi_ttt, hp, hc);
314 {
printf(
"%.6le %.6le %.6le %.6le %.6le %.6le\n", r,
pr, phi,
tau,
r2,
a2);
317 double r_tt, r_ttt, phi_tt, phi_ttt;
321 printf(
"%.6le %.6le %.6le %.6le %.6le\n", sqrt(
a2), r_tt, r_ttt, phi_tt, phi_ttt);
322 h_quad(::r,
r_dot, r_tt, r_ttt, ::phi,
w, phi_tt, phi_ttt, hp, hc);
323 printf(
"%.6le %.6le\n", hp, hc);
wavearray< double > t(hp.size())
void geo_higher_deriv(double E, double L, double a, double &r_dotdot, double &r_dotdotdot, double &wdot, double &wdotdot)
double min(double x, double y)
printf("total live time: non-zero lags = %10.1f \n", liveTot)
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
geodesic(double r0, double phi0, double E0, double L0, int dir, double Jspin, double mu=0, double dt=1)
bool integrate(double D_t, int &N, double rmin=0, bool stop_alr=false)
void reduced_quad_dot(double r, double r_t, double r_tt, double r_ttt, double phi, double phi_t, double phi_tt, double phi_ttt, double(&I_tt)[3][3], double(&I_ttt)[3][3])
void printfutil(double *y)
double fabs(const Complex &x)
void h_quad(double r, double r_t, double r_tt, double r_ttt, double phi, double phi_t, double phi_tt, double phi_ttt, double &hplus, double &hcross)
void f_util(double *y, double Jspin, double mu)
int f(realtype t, N_Vector y, N_Vector ydot, void *ud)
void quad_loss(double r, double r_t, double r_tt, double r_ttt, double phi, double phi_t, double phi_tt, double phi_ttt, double &E_dot, double &Lz_dot)