From 2c05355eaab505a003b94aa4891a193b93e41183 Mon Sep 17 00:00:00 2001 From: Yohai Meiron Date: Fri, 13 Mar 2020 16:58:41 -0400 Subject: [PATCH] Removed more ifdef blocks --- act_def_linklist.c | 376 ------------------- phigrape.cpp | 889 +++------------------------------------------ pn_bh.c | 765 -------------------------------------- star_destr.c | 275 -------------- star_destr_ext.c | 286 --------------- 5 files changed, 56 insertions(+), 2535 deletions(-) delete mode 100644 act_def_linklist.c delete mode 100644 pn_bh.c delete mode 100644 star_destr.c delete mode 100644 star_destr_ext.c diff --git a/act_def_linklist.c b/act_def_linklist.c deleted file mode 100644 index b879bc3..0000000 --- a/act_def_linklist.c +++ /dev/null @@ -1,376 +0,0 @@ -/************************************************************** - File : act_def_linklist.c - Func. : provide linear linking list functions - : for active particle def. - CODED BY : Zhong Shiyan - START : 2014-03-28, 12:30 -**************************************************************/ - -//***********************************************************// -/* Definition of T/P-node */ -//***********************************************************// -typedef struct PNODE -{ - int Pid; // Particle's real ID - struct PNODE *NextPNODE; -} PNODE; - -typedef struct TNODE -{ - double t_node; // t_node = t + dt - int n_node; // number of P-nodes under this node - - struct PNODE *PartList, *PartListEnd; - struct TNODE *NextTNODE; -} TNODE; - -struct TNODE *CurrT=NULL; - -//***********************************************************// -/* Operations on T/P-node */ -//***********************************************************// -struct TNODE *CreateTNODE( double t ){ - - struct TNODE *ptr; - - ptr = (TNODE*)malloc(sizeof(*ptr)); - if( ptr == NULL ){ - printf("Fail to create a node."); - exit(-1); - } - - ptr->t_node = t; - ptr->n_node = 0; - - ptr->NextTNODE = NULL; - - ptr->PartList = NULL; - ptr->PartListEnd = NULL; - - return ptr; -} -//***********************************************************// -struct PNODE *DeletePNODE( struct PNODE *ptr ){ - struct PNODE *next; - - next = ptr->NextPNODE; - free(ptr); - - return next; -} -//***********************************************************// -void DeleteTNODE( struct TNODE *Tptr ){ - - struct PNODE *Pptr; - - Pptr = Tptr->PartList; - - while( Pptr != NULL ){ - Pptr = DeletePNODE( Pptr ); - } - - free(Tptr); - return; -} -//***********************************************************// -struct PNODE *CreatePNODE( int id ){ - struct PNODE *ptr; - - ptr = (PNODE*)malloc(sizeof(*ptr)); - if( ptr == NULL ){ - printf("Fail to create a P-node."); - return NULL; - } - - ptr->Pid = id; - ptr->NextPNODE = NULL; - - return ptr; - -} -//***********************************************************// -void InsertTNODE( struct TNODE *front, - struct TNODE *rear, - struct TNODE *ptr ) -{ - front->NextTNODE = ptr; - ptr->NextTNODE = rear; - -} -//***********************************************************// - - - - -//***********************************************************// -/* Functions about act_def */ -//***********************************************************// -/************************************************************** -This function only called 1 time, at beginning of simulation. -After all particle's dt are computed -**************************************************************/ - -void CreateLinkList(){ - - struct TNODE *head,*tail,*Tp1,*Tp2,*newTNODE,*Tptr; - struct PNODE *Pptr,*Ptail; - - int i, iii; - double ttmp, t1, t2; - - head = CreateTNODE( -1.0 ); // will be discarded at the end of this routine - tail = CreateTNODE( 2.0*t_end ); - - head->NextTNODE = tail; - tail->NextTNODE = NULL; - - // building link list - -for(i=0;iNextTNODE; - t1 = Tp1->t_node; - t2 = Tp2->t_node; - - Pptr = CreatePNODE( iii ); - - while( Tp1->NextTNODE != NULL ){ - if( ttmp == t1 ){ // if T-node exist, add a P-node - - if( Tp1->PartListEnd == NULL ){ // This is the first P-node under current T-node - Tp1->PartList = Pptr; - Tp1->PartListEnd = Pptr; - Tp1->n_node = Tp1->n_node + 1; - } - else{ // There are already many P-nodes under this T-node - Ptail = Tp1->PartListEnd; - Ptail->NextPNODE = Pptr; - Tp1->PartListEnd = Pptr; - Tp1->n_node = Tp1->n_node + 1; - } - break; // jump out of this "while" loop - } - - if( ttmp > t1 && ttmp < t2 ){ // Create a new T-node and insert between *Tp1 and *Tp2, then add P-node to it - - newTNODE = CreateTNODE(ttmp); - InsertTNODE( Tp1, Tp2, newTNODE); - - newTNODE->n_node = 1; - newTNODE->PartList = Pptr; - newTNODE->PartListEnd = Pptr; - break; // jump out of this "while" loop - } - - // move to next T-node - Tp1 = Tp1->NextTNODE; - t1 = Tp1->t_node; - - if(Tp2->NextTNODE != NULL){ - Tp2 = Tp2->NextTNODE; - t2 = Tp2->t_node; - } - else{ break; } - - }// while( Tp1->NextTNODE != NULL ) - -}//for(i=0;iNextTNODE; - free(head); - -} - -//End of CreateLinkList() -/**************************************************************/ - - -/************************************************************** -This Function is used to modify the link list, -after get new dt for active particles. Then point *CurrT to next -T-node. -**************************************************************/ - -void ModifyLinkList(){ - - struct TNODE *Tp1,*Tp2,*newTNODE; - struct PNODE *Pptr,*Ptail; - - int i, iii; - double ttmp, t1, t2; - - for(i=0;iNextTNODE; - t1 = Tp1->t_node; - t2 = Tp2->t_node; - - Pptr = CreatePNODE( iii ); - - while( Tp1->NextTNODE != NULL ){ - if( ttmp == t1 ){ // if T-node exist, add a P-node - - if( Tp1->PartListEnd == NULL ){ // This is the first P-node under current T-node - Tp1->PartList = Pptr; - Tp1->PartListEnd = Pptr; - Tp1->n_node = Tp1->n_node + 1; - } - else{ // There are already many P-nodes under this T-node - Ptail = Tp1->PartListEnd; - Ptail->NextPNODE = Pptr; - Tp1->PartListEnd = Pptr; - Tp1->n_node = Tp1->n_node + 1; - } - break; // jump out of this "while" loop - } - - if( ttmp > t1 && ttmp < t2 ){ // Create a new T-node and insert between *Tp1 and *Tp2, then add P-node to it - - newTNODE = CreateTNODE(ttmp); - InsertTNODE( Tp1, Tp2, newTNODE); - - newTNODE->n_node = 1; - newTNODE->PartList = Pptr; - newTNODE->PartListEnd = Pptr; - - break; // jump out of this "while" loop - } - - // move to next T-node - Tp1 = Tp1->NextTNODE; - t1 = Tp1->t_node; - - if(Tp2->NextTNODE != NULL){ - Tp2 = Tp2->NextTNODE; - t2 = Tp2->t_node; - } - else{ break; } - - }//while( Tp1->NextTNODE != NULL ) - - }//for(i=0;iNextTNODE; - DeleteTNODE( Tp1 ); -} - -//End of ModifyLinkList() -/***************************************************************/ - - -/***************************************************************/ -/* -void i_swap(int *a, int *b) -{ -register int tmp; -tmp = *a; *a = *b; *b = tmp; -} -*/ -/***************************************************************/ -/***************************************************************/ -void ind_act_sort(int l, int r) -{ - -int i, j, cikl, tmp; - -i = l; j = r; -tmp = ind_act[(l+r)/2]; - -cikl = 1; - -while(cikl) - { - while (ind_act[i]PartList; - n_act = 0; - - min_t = CurrT->t_node; // IMPORTANT !! - - flag = 0; - - while(Pptr != NULL) - { - iii = Pptr->Pid; - Pptr = Pptr->NextPNODE; - if( m[iii] != 0.0 ) // Do not put zero mass part. in active plist - { - ind_act[i] = iii; -// if(ind_act[i]==N-1) flag=1; - i++; n_act++; - } - } - - if( n_act > 2 ) ind_act_sort( 0, n_act-1 ); - -// printf("last pid: %06d\n", ind_act_ll[n_act-1]); -// if( flag == 0 ) printf("Warning: BH not in the ind_act array!\n"); -// if( ind_act[n_act-1]!= N-1) printf("Warning: BH not in the last of ind_act array!\n"); - - }// End of get_act_plist() -/**************************************************************/ - - -#ifdef DEBUG_extra -/************************************************************** -Check link list -**************************************************************/ -void check_linklist(int Tstep) - { - FILE *listf; - struct TNODE *Tptr; - - listf = fopen("Check_linklist.dat","a"); - Tptr = CurrT; - fprintf(listf,"Timesteps = %04d\n", Tstep); - while(Tptr != NULL) - { - fprintf(listf,"% 8E %04d\n", Tptr->t_node, Tptr->n_node); - Tptr = Tptr->NextTNODE; - } - - fprintf(listf,"========================\n\n"); - - fclose(listf); - } -/**************************************************************/ -#endif - diff --git a/phigrape.cpp b/phigrape.cpp index d258c08..70804df 100644 --- a/phigrape.cpp +++ b/phigrape.cpp @@ -137,7 +137,6 @@ Last redaction : 2019.04.16 12:55 #error "Contradicting preprocessor flags!" #endif -/****************************************************************************/ #include #include #include @@ -231,8 +230,6 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, E_tot_corr, E_tot_corr_0, DE_tot_corr, E_tot_corr_sd, E_tot_corr_sd_0, DE_tot_corr_sd, E_corr = 0.0, E_sd = 0.0, t_diss_on = 0.125, - e_kin_corr = 0.0, e_pot_corr = 0.0, e_pot_BH_corr = 0.0, - e_summ_corr = 0.0, mcm, rcm_mod, vcm_mod, rcm_sum=0.0, vcm_sum=0.0, eps=0.0, eps2, @@ -252,8 +249,6 @@ double dt_disk, dt_contr, t_disk=0.0, t_contr=0.0, tmp_a_bh_pn0, tmp_a_bh_pn1, tmp_a_bh_pn2, tmp_a_bh_pn2_5, tmp_a_bh_pn3, tmp_a_bh_pn3_5, tmp_a_bh_spin, tmp_adot_bh_pn0, tmp_adot_bh_pn1, tmp_adot_bh_pn2, tmp_adot_bh_pn2_5, tmp_adot_bh_pn3, tmp_adot_bh_pn3_5, tmp_adot_bh_spin; -double R[3], V[3], L[3], mju, dr2, dr, dv2, dv, dl2, dl, pot_eff; - char processor_name[MPI_MAX_PROCESSOR_NAME], inp_fname[30], out_fname[30], dbg_fname[30]; @@ -335,10 +330,6 @@ double eps_BH=0.0; double m_bulge, a_bulge, b_bulge, m_disk, a_disk, b_disk, m_halo, a_halo, b_halo, -// x_ext, y_ext, z_ext, -// vx_ext, vy_ext, vz_ext, -// x_ij, y_ij, z_ij, -// vx_ij, vy_ij, vz_ij, rv_ij, x2_ij, y2_ij, z2_ij, r_tmp, r2_tmp, z_tmp, z2_tmp; #endif @@ -359,12 +350,6 @@ double r2, rv_ij, m_bh, b_bh, eps_bh; #endif -#ifdef EXTPOT_SC -int M_SC=M_SC_DIM; -double r_sc[M_SC_DIM], m_sc[M_SC_DIM], p_sc[M_SC_DIM]; -double M_R, pot_out_R; -#endif - double pot_ext[N_MAX], pot_act_ext[N_MAX]; // for all EXTPOT #endif @@ -383,9 +368,6 @@ double x_bbhc[3], v_bbhc[3], DR2, tmp_r2; double DV2, EB, SEMI_a, SEMI_a2; #endif -/****************************************************************************/ - -/****************************************************************************/ #ifdef ADD_N_BH double x_bh1[3], x_bh2[3], v_bh1[3], v_bh2[3]; @@ -432,9 +414,7 @@ return - 0 if everything OK */ #endif -/****************************************************************************/ -/****************************************************************************/ #ifdef ADD_PN_BH double C_NB = 477.12; @@ -447,11 +427,7 @@ double a_pn1[7][3], adot_pn1[7][3], double s_bh1[3] = {0.0, 0.0, 1.0}; double s_bh2[3] = {0.0, 0.0, 1.0}; -#ifdef ADD_SPIN_BH // eto rabotajet vsegda !!! #include "pn_bh_spin.c" -#else -#include "pn_bh.c" // eto staraja versija, bolshe ne rabotajet !!! -#endif /* int calc_force_pn_BH(double m1, double xx1[], double vv1[], double ss1[], @@ -495,9 +471,6 @@ return - 0 if everything OK #endif -/****************************************************************************/ - -/****************************************************************************/ /* RAND_MAX = 2147483647 */ /* my_rand : 0.0 - 1.0 */ /* my_rand2 : -1.0 - 1.0 */ @@ -511,19 +484,6 @@ double my_rand2(void) { return (double)(2.0)*((rand() - RAND_MAX/2)/(double)RAND_MAX); } -/****************************************************************************/ - -#ifdef STARDESTR -#include "star_destr.c" -#endif - -#ifdef STARDESTR_EXT -#include "star_destr_ext.c" -#endif - -#ifdef ACT_DEF_LL -#include "act_def_linklist.c" -#endif #ifdef ETICS double t_exp, dt_exp=ETICS_DTSCF; // t_exp is just the initial value @@ -532,7 +492,6 @@ int grapite_cep_index; #endif #endif -/****************************************************************************/ void get_CPU_time(double *time_real, double *time_user, double *time_syst) { struct rusage xxx; @@ -558,9 +517,7 @@ gettimeofday(&tv, NULL); *time_user = *time_real; } -/****************************************************************************/ -/****************************************************************************/ void read_data() { inp = fopen(inp_fname,"r"); @@ -574,34 +531,25 @@ for (i=0; i r) ) - { - r1 = r_sc[ii]; r2 = r_sc[ii+1]; - m1 = m_sc[ii]; m2 = m_sc[ii+1]; - p1 = p_sc[ii]; p2 = m_sc[ii+1]; - break; - } - } - -m_tmp = m1 + (r-r1)*(m2-m1)/(r2-r1); -p_tmp = p1 + (r-r1)*(p2-p1)/(r2-r1); - -// if (myRank == rootRank) -// { -// printf("\t \t \t %.6E %.6E \n", r1, r2); -// printf("\t \t \t %.6E %.6E \n", m1, m2); -// printf("\t \t \t %.6E %.6E \n", p1, p2); -// printf("\t \t \t %06d %.6E %.6E %.6E \n", ii, r, m_tmp, p_tmp); -// fflush(stdout); -// } /* if (myRank == rootRank) */ - -*M_R = m_tmp; -*pot_out_R = p_tmp; - -//exit(-1); - -} -#endif -/****************************************************************************/ - -/****************************************************************************/ void calc_ext_grav_zero() { @@ -1492,171 +1292,73 @@ for (i=0; i 100.0) jerk_i[ii][0] = 100.0; - if (ABS(jerk_i[ii][1]) > 100.0) jerk_i[ii][1] = 100.0; - if (ABS(jerk_i[ii][2]) > 100.0) jerk_i[ii][2] = 100.0; - } - - g6calc_firsthalf(clusterid, n_loc, nn, index_i, x_i, v_i , a_i, jerk_i, p_i, eps2, h2_i); - g6calc_lasthalf(clusterid, n_loc, nn, index_i, x_i, v_i, eps2, h2_i, a_i, jerk_i, p_i); - - g6_calls++; -*/ - - for (ii=0; ii R_LIMITS_ALL) && (tmp_rv > 0.0) ) { - for (k=0;k<3;k++) { - v[i][k] *= 1.0E-01; - } /* k */ - } /* if */ - } /* i */ -#endif - -#ifdef CMCORR - /* possible CM correction of the initial datafile... */ - if ( (diskstep == 0) && (time_cur == 0.0) ) cm_corr(); -#endif printf("\n"); printf("Begin the calculation of phi-GRAPE program on %03d processors\n", n_proc); printf("\n"); @@ -2327,12 +1992,6 @@ int main(int argc, char *argv[]) printf("v_halo = %.6E r_halo = %.6E \n", v_halo, r_halo); #endif -#ifdef EXTPOT_SC - read_SC_mass(); - - printf("EXTPOT_SC:# of points %06d \t M_SC_TOT = %.6E \n", M_SC_DIM-1, m_sc[M_SC_DIM-1]); -#endif - printf("\n"); #endif @@ -2377,104 +2036,12 @@ int main(int argc, char *argv[]) #endif // ADD_BH2 -#ifdef STARDESTR - -#ifdef ADD_BH2 - out = fopen("mass-bh.dat","w"); - m_bh1 = m[0]; - m_bh2 = m[1]; - num_bh1 = 0; - num_bh2 = 0; - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - - out = fopen("mass-bh.con","w"); - m_bh1 = m[0]; - m_bh2 = m[1]; - num_bh1 = 0; - num_bh2 = 0; - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); -#else - -#ifdef ADD_BH1 - - out = fopen("mass-bh.dat","w"); - m_bh1 = m[0]; - num_bh1 = 0; - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - - out = fopen("mass-bh.con","w"); - m_bh1 = m[0]; - num_bh1 = 0; - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#endif // STARDESTR - -#ifdef STARDESTR_EXT - num_bh = 0; - out = fopen("mass-bh.dat","w"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - - num_bh = 0; - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - -#ifdef BH_OUT_NB_EXT - out = fopen("bh_nb.dat","w"); - fclose(out); -#endif - -#endif // STARDESTR_EXT - #ifdef BBH_INF out = fopen("bbh.inf","w"); fclose(out); #endif } else { // if (diskstep == 0) -#ifdef STARDESTR - -#ifdef ADD_BH2 - inp = fopen("mass-bh.con","r"); - fscanf(inp,"%lE %lE %d %lE %d", &tmp, &m_bh1, &num_bh1, &m_bh2, &num_bh2); - fclose(inp); - printf("%.8E \t %.8E %06d \t %.8E %06d \n", tmp, m_bh1, num_bh1, m_bh2, num_bh2); - fflush(stdout); -#else - -#ifdef ADD_BH1 - inp = fopen("mass-bh.con","r"); - fscanf(inp,"%lE %lE %d", &tmp, &m_bh1, &num_bh1); - fclose(inp); - printf("%.8E \t %.8E \t %06d \n", tmp, m_bh1, num_bh1); - fflush(stdout); -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#endif // STARDESTR - -#ifdef STARDESTR_EXT - inp = fopen("mass-bh.con","r"); - fscanf(inp,"%lE %lE %d", &tmp, &m_bh, &num_bh); - fclose(inp); - printf("%.8E \t %.8E \t %06d \n", tmp, m_bh, num_bh); - fflush(stdout); -#endif // STARDESTR_EXT - } // if (diskstep == 0) get_CPU_time(&CPU_time_real0, &CPU_time_user0, &CPU_time_syst0); @@ -2539,18 +2106,6 @@ int main(int argc, char *argv[]) #endif // EXTPOT -#ifdef STARDESTR - MPI_Bcast(&m_bh1, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&num_bh1, 1, MPI_INT, rootRank, MPI_COMM_WORLD); - - MPI_Bcast(&m_bh2, 1, MPI_DOUBLE, rootRank, MPI_COMM_WORLD); - MPI_Bcast(&num_bh2, 1, MPI_INT, rootRank, MPI_COMM_WORLD); -#endif // STARDESTR - -#ifdef STARDESTR_EXT - MPI_Bcast(&num_bh, 1, MPI_INT, rootRank, MPI_COMM_WORLD); -#endif // STARDESTR_EXT - /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); @@ -2867,14 +2422,6 @@ int main(int argc, char *argv[]) if (dt_tmp < dt_min) dt_tmp = dt_min; if (dt_tmp > dt_max) dt_tmp = dt_max; -#ifdef STARDESTR - if (m[i] == 0.0) dt_tmp = dt_max; -#endif - -#ifdef STARDESTR_EXT - if (m[i] == 0.0) dt_tmp = dt_max; -#endif - dt[i] = dt_tmp; #ifdef DT_MIN_WARNING @@ -2919,10 +2466,6 @@ int main(int argc, char *argv[]) #endif // ADD_BH2 -#ifdef ACT_DEF_LL - CreateLinkList(); -#endif - /* Wait to all processors to finish his works... */ MPI_Barrier(MPI_COMM_WORLD); @@ -3008,21 +2551,6 @@ int main(int argc, char *argv[]) /* Define the minimal time and the active particles on all the nodes (exclude the ZERO masses!!!) */ -#ifdef ACT_DEF_LL - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); -#endif - - get_act_plist(); - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real, &CPU_tmp_user, &CPU_tmp_syst); - DT_ACT_DEF1 += (CPU_tmp_user - CPU_tmp_user0); -#endif - -#else - #ifdef TIMING get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); #endif @@ -3117,8 +2645,6 @@ int main(int argc, char *argv[]) DT_ACT_DEF2 += (CPU_tmp_user - CPU_tmp_user0); #endif -#endif // ACT_DEF_LL - #ifdef TIMING get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); #endif @@ -3344,14 +2870,6 @@ int main(int argc, char *argv[]) #endif // ADD_BH2 -#ifdef STARDESTR - if (m_act[i] == 0.0) dt_new = dt_max; -#endif - -#ifdef STARDESTR_EXT - if (m_act[i] == 0.0) dt_new = dt_max; -#endif - if (dt_new < dt_min) dt_tmp = dt_min; if ( (dt_new < dt_tmp) && (dt_new > dt_min) ) { @@ -3519,126 +3037,6 @@ int main(int argc, char *argv[]) DT_ACT_CORR += (CPU_tmp_user - CPU_tmp_user0); #endif - /* STARDESTR routine on all the nodes */ - -#ifdef STARDESTR - -#ifdef TIMING - get_CPU_time(&CPU_tmp_real0, &CPU_tmp_user0, &CPU_tmp_syst0); -#endif - -#ifdef ADD_BH2 - - star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act, - N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1); - - m_act[i_bh1] = m_bh1; - - star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act, - N, ind, m, x, v, pot, a, adot, t, &m_bh2, &num_bh2, i_bh2); - - m_act[i_bh2] = m_bh2; - -#else - -#ifdef ADD_BH1 - - star_destr(min_t, n_act, ind_act, m_act, x_act, v_act, pot_act, a_act, adot_act, t_act, dt_act, - N, ind, m, x, v, pot, a, adot, t, &m_bh1, &num_bh1, i_bh1); - - m_act[i_bh1] = m_bh1; - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#ifdef LIMITS - for (i=0; i= t_contr) { if (myRank == rootRank) { -#ifdef STARDESTR - -#ifdef ADD_BH2 - - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E %06d \t %.8E %06d \n", - time_cur, m_bh1, num_bh1, m_bh2, num_bh2); - fclose(out); - -#else - -#ifdef ADD_BH1 - - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E %06d \n", - time_cur, m_bh1, num_bh1); - fclose(out); - -#endif // ADD_BH1 - -#endif // ADD_BH2 - -#endif // STARDESTR - -#ifdef STARDESTR_EXT - - out = fopen("phi-GRAPE.ext","w"); - fprintf(out,"%.8E \t %.8E \n", m_bh, b_bh); - fclose(out); - - out = fopen("mass-bh.dat","a"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - - out = fopen("mass-bh.con","w"); - fprintf(out,"%.8E \t %.8E \t %06d \n", time_cur, m_bh, num_bh); - fclose(out); - -#endif // STARDESTR_EXT energy_contr(); -#ifdef CMCORR111 - for (i=0; i R_LIMITS_ALL) && (tmp_rv > 0.0) ) { - for (k=0;k<3;k++) { - v[i][k] *= 1.0E-01; - } /* k */ - } /* if */ - } /* i */ - - for (j=0; j 900.0 ) - { - - m[i] = 0.0; - - pot[i] = 0.0; - -#ifdef EXTPOT - pot_ext[i] = 0.0; -#endif - - for (k=0;k<3;k++) - { - x[i][k] = 1000.0 + 1.0*my_rand2(); - v[i][k] = 1.0E-06*my_rand2(); - a[i][k] = 1.0E-06*my_rand2(); - adot[i][k] = 1.0E-06*my_rand2(); - } /* k */ - - t[i] = 2.0*t_end; - - dt[i] = 0.125; - - } /* if */ - - } /* i */ - - for (j=0; j= t_disk) */ -#ifdef CPU_TIMELIMIT - if (myRank == rootRank) { - get_CPU_time(&CPU_time_real, &CPU_time_user, &CPU_time_syst); - tmp_cpu = CPU_time_real-CPU_time_real0; - } /* if (myRank == rootRank) */ -#endif - -#ifdef ACT_DEF_LL - ModifyLinkList(); -#endif } /* while (time_cur < t_end) */ /* close the local GRAPE's */ diff --git a/pn_bh.c b/pn_bh.c deleted file mode 100644 index 54822f3..0000000 --- a/pn_bh.c +++ /dev/null @@ -1,765 +0,0 @@ -/***************************************************************************/ -/* - Coded by : Peter Berczik (on the base of Gabor Kupi original PN code) - Version number : 2.0 SPIN - Last redaction : 2012.V.07. 11:16 -*/ - -int calc_force_pn_BH(double m1, double xx1[], double vv1[], double spin1[], - double m2, double xx2[], double vv2[], double spin2[], - double CCC_NB, double dt_bh, - int usedOrNot[], - double a_pn1[][3], double adot_pn1[][3], - double a_pn2[][3], double adot_pn2[][3]) -{ - -/* - INPUT - -m1 - mass of the 1 BH -xx1[0,1,2] - coordinate of the 1 BH -vv1[0,1,2] - velocity of the 1 BH -spin1[0,1,2] - normalized spin of the 1 BH - -m2 - mass of the 2 BH -xx2[0,1,2] - coordinate of the 2 BH -vv2[0,1,2] - velocity of the 2 BH -spin2[0,1,2] - normalized spin of the 2 BH - -CCC_NB - Speed of light "c" in internal units -dt_BH - timestep of the BH's, needed for the SPIN integration - -usedOrNot[PN0, PN1, PN2, PN2.5, PN3, PN3.5, SPIN] - different PN term usage: PN1, PN2, PN2.5, PN3, PN3.5, SPIN - 0 1 2 3 4 5 6 - - OUTPUT - -a_pn1 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH -adot_pn1[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 1 BH - -a_pn2 [0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH -adot_pn2[0 - PN0; 1 - PN1; 2 - PN2; 3 - PN2.5, 4 - PN3, 5 - PN3.5, 6 - SPIN] [3] for the 2 BH - -return - 0 if everything OK - - 505 if BH's separation < 4 x (RSwarch1 + RSwarch2) -*/ - - - -int j, k; -double PI2 = 9.86960440108935; - -double c_1, c_2, c_4, c_5, c_6, c_7, RS_DIST; - -double M, eta, r, r2, r3, MOR; -double V1_V22,VWHOLE, RP, RPP, VA; -double N[3], x[3], v[3], A[3]; -double A1, B1, A2, B2, A2_5, B2_5, AK2, BK2, AK4, BK4, AK5, BK5; -double A1D, A2D, A2_5D, B1D, B2D, B2_5D, ADK2, BDK2, ADK4, BDK4, ADK5, BDK5; - -double A3, B3, A3_5, B3_5, AK6, BK6, AK7, BK7; -double A3D, A3_5D, B3D, B3_5D, ADK6, BDK6, ADK7, BDK7; - -int Van_Spin=0; -int Van_QM=0; - -double DM, S1[3], SPIN[3][2], S2[3], KSS[3], KSSIG[3], XS[3], XA[3], NCV[3], NCS[3], NCSIG[3], - VCS[3], VCSIG[3], SDNCV, SIGDNCV, NDV, XS2, XA2, NXA, NXS, VDS, VDSIG, NDS, NDSIG, - C1_5[3], C2[3], C2_5[3]; -double LABS, LU[3], S1DLU, S2DLU, SU1[3], SV1[3], SS1[3], SS2[3], SU2[3], SV2[3]; -double AT[3], NDOT[3], NVDOT, NDOTCV[3], NCA[3]; -double SS1aux[3],SS2aux[3],SU[3],SV[3],XAD[3],XSD[3]; -double NDOTCS[3], NCSU[3], NDOTCSIG[3], NCSV[3], ACS[3], VCSU[3], ACSIG[3], VCSV[3], SNVDOT, - SIGNVDOT, NSDOT, NSIGDOT, VSDOT, VSIGDOT, NXSDOT, NXADOT; -double C1_5D[3],C2D[3], C2_5D[3]; -double ADK, BDK, AD[3], KSAK, KSBK; -double nu, Spin1Abs2, Spin2Abs2, rS1, rS2, S1Dir[3], S2Dir[3], QM[3]; -double Spin1Abs, Spin2Abs, QMAux2_1[3], QMAux2_2[3], QMAux1[3] , QMD[3], SPINPrev[3][2], SpinPrev2_1, - SpinPrev2_2, SPSPP1, SPSPP2, Spin1AbsNew2, Spin2AbsNew2, Spin1AbsNew, Spin2AbsNew, S1DirNew[3], - S2DirNew[3], rS1p, rS2p, S1p[3], S2p[3], Np[3]; - - -Van_Spin = usedOrNot[6]; // Van vagy nincs SPIN szamolas... - - -for(k=0;k<3;k++) - { - SPIN[k][0] = spin1[k]; - SPIN[k][1] = spin2[k]; - } - - -for(j=0;j<7;j++) - { -for(k=0;k<3;k++) - { - a_pn1[j][k] = 0.0; adot_pn1[j][k] = 0.0; - a_pn2[j][k] = 0.0; adot_pn2[j][k] = 0.0; - } - } - - -// Speed of light "c" and its powers - -c_1 = CCC_NB; - -c_2 = SQR(c_1); -c_4 = SQR(c_2); -c_5 = c_4*c_1; -c_6 = c_5*c_1; -c_7 = c_6*c_1; - - -// Mass parameters - -M = m1+m2; -eta = m1*m2/(M*M); -nu = m1/m2; - -for(k=0;k<3;k++) - { - x[k] = xx1[k] - xx2[k]; - v[k] = vv1[k] - vv2[k]; - } - - -r2 = SQR(x[0]) + SQR(x[1]) + SQR(x[2]); -r = sqrt(r2); -r3 = r2*r; - -MOR = M/r; -V1_V22 = v[0]*v[0]+v[1]*v[1]+v[2]*v[2]; -VWHOLE = sqrt(V1_V22); -RP = (x[0]*v[0]+x[1]*v[1]+x[2]*v[2])/r; - - -// Newton accelerations - -for(k=0;k<3;k++) N[k] = x[k]/r; - -// PN accelerations - -AK2 = 0.0; BK2 = 0.0; -AK4 = 0.0; BK4 = 0.0; -AK5 = 0.0; BK5 = 0.0; -AK6 = 0.0; BK6 = 0.0; -AK7 = 0.0; BK7 = 0.0; - -for(k=0;k<3;k++) - { - C1_5[k] = 0.0; - C2[k] = 0.0; - C2_5[k] = 0.0; - QM[k] = 0.0; - } - - -if(usedOrNot[1] == 1) // PN1 ~1/c^2 - { - A1 = 2.0*(2.0+eta)*MOR-(1.0+3.0*eta)*V1_V22 +1.5*eta*RP*RP; - B1 = 2.0*(2.0-eta)*RP; - - AK2 = A1/c_2; - BK2 = B1/c_2; - } - -if(usedOrNot[2] == 1) // PN2 ~1/c^4 - { - A2 = -0.75*(12.0+29.0*eta)*MOR*MOR-eta*(3.0-4.0*eta)*V1_V22*V1_V22-1.875*eta*(1.0-3.0*eta)*RP*RP*RP*RP+0.5*eta*(13.0-4.0*eta)*MOR*V1_V22+(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RP+1.5*eta*(3.0-4.0*eta)*V1_V22*RP*RP; - B2 = -0.5*RP*((4.0+41.0*eta+8.0*eta*eta)*MOR-eta*(15.0+4.0*eta)*V1_V22+3.0*eta*(3.0+2.0*eta)*RP*RP); - - AK4 = A2/c_4; - BK4 = B2/c_4; - } - -if(usedOrNot[3] == 1) // PN2.5 ~1/c^5 - { - A2_5 = 1.6*eta*MOR*RP*(17.0*MOR/3.0+3.0*V1_V22); - B2_5 = -1.6*eta*MOR*(3.0*MOR+V1_V22); - - AK5 = A2_5/c_5; - BK5 = B2_5/c_5; - } - -if(usedOrNot[4] == 1) // PN3 ~1/c^6 - { - A3 = MOR*MOR*MOR*(16.0+(1399.0/12.0-41.0*PI2/16.0)*eta+ - 71.0*eta*eta/2.0)+eta*(20827.0/840.0+123.0*PI2/64.0-eta*eta) - *MOR*MOR*V1_V22-(1.0+(22717.0/168.0+615.0*PI2/64.0)*eta+ - 11.0*eta*eta/8.0-7.0*eta*eta*eta)*MOR*MOR*RP*RP- - 0.25*eta*(11.0-49.0*eta+52.0*eta*eta)*V1_V22*V1_V22*V1_V22+ - 35.0*eta*(1.0-5.0*eta+5.0*eta*eta)*RP*RP*RP*RP*RP*RP/16.0- - 0.25*eta*(75.0+32.0*eta-40.0*eta*eta)*MOR*V1_V22*V1_V22- - 0.5*eta*(158.0-69.0*eta-60.0*eta*eta)*MOR*RP*RP*RP*RP+ - eta*(121.0-16.0*eta-20.0*eta*eta)*MOR*V1_V22*RP*RP+ - 3.0*eta*(20.0-79.0*eta+60.0*eta*eta)*V1_V22*V1_V22*RP*RP/8.0- - 15.0*eta*(4.0-18.0*eta+17.0*eta*eta)*V1_V22*RP*RP*RP*RP/8.0; - - B3 = RP*((4.0+(5849.0/840.0+123.0*PI2/32.0)*eta-25.0*eta*eta- - 8.0*eta*eta*eta)*MOR*MOR+eta*(65.0-152.0*eta-48.0*eta*eta)* - V1_V22*V1_V22/8.0+15.0*eta*(3.0-8.0*eta-2.0*eta*eta)*RP*RP*RP*RP/8.0+ - eta*(15.0+27.0*eta+10.0*eta*eta)*MOR*V1_V22-eta*(329.0+177.0*eta+ - 108.0*eta*eta)*MOR*RP*RP/6.0- - 3.0*eta*(16.0-37.0*eta-16.0*eta*eta)*V1_V22*RP*RP/4.0); - - AK6 = A3/c_6; - BK6 = B3/c_6; - } - -if(usedOrNot[5] == 1) // PN3.5 ~1/c^7 - { - A3_5 = MOR*eta*(V1_V22*V1_V22*(-366.0/35.0-12.0*eta)+V1_V22*RP*RP*(114.0+12.0*eta)-112.0*RP*RP*RP*RP+MOR*(V1_V22*(-692.0/35.0+724.0*eta/15.0)+RP*RP*(-294.0/5.0-376.0*eta/5.0)+MOR*(-3956.0/35.0-184.0*eta/5.0))); - B3_5 = 8.0*eta*MOR*((1325.0+546.0*eta)*MOR*MOR/42.0+(313.0+42.0*eta)*V1_V22*V1_V22/28.0+75.0*RP*RP*RP*RP-(205.0+777.0*eta)*MOR*V1_V22/42.0+(205.0+424.0*eta)*MOR*RP*RP/12.0-3.0*(113.0+2.0*eta)*V1_V22*RP*RP/4.0)/5.0; - - AK7 = A3_5/c_7; - BK7 = B3_5/c_7; - } - - -// Spin accelerations - - if(Van_Spin==1) - { - DM = m1 - m2; - - Spin1Abs2 = 0.0; - Spin2Abs2 = 0.0; - rS1 = 0.0; - rS2 = 0.0; - - for(k=0;k<3;k++) - { - Spin1Abs2 += SPIN[k][0]*SPIN[k][0]; // normalizalt spin - Spin2Abs2 += SPIN[k][1]*SPIN[k][1]; - rS1 += N[k]*S1Dir[k]; - rS2 += N[k]*S2Dir[k]; - - S1[k] = SPIN[k][0]*m1*m1/c_1; // fizikai spin - S2[k] = SPIN[k][1]*m2*m2/c_1; - KSS[k] = S1[k]+S2[k]; - KSSIG[k] = M*(S2[k]/m2-S1[k]/m1); - XS[k] = 0.5*(SPIN[k][0]+SPIN[k][1]); - XA[k] = 0.5*(SPIN[k][0]-SPIN[k][1]); - } - - Spin1Abs = sqrt(Spin1Abs2); - Spin2Abs = sqrt(Spin2Abs2); - - for(k=0;k<3;k++) - { - S1Dir[k] = SPIN[k][0]/Spin1Abs; - S2Dir[k] = SPIN[k][1]/Spin2Abs; - } - - - //NCV crossproduct of N[k] and relative v = N[k]Xv[j] - NCV[0] = N[1]*v[2] - N[2]*v[1]; - NCV[1] = N[2]*v[0] - N[0]*v[2]; - NCV[2] = N[0]*v[1] - N[1]*v[0]; - - //NCS crossproduct of N[k] and KSS = N[k]XKSS - NCS[0] = N[1]*KSS[2] - N[2]*KSS[1]; - NCS[1] = N[2]*KSS[0] - N[0]*KSS[2]; - NCS[2] = N[0]*KSS[1] - N[1]*KSS[0]; - - //NCSIG crossproduct of N[k] and KSSIG = N[k]XKSSIG - NCSIG[0] = N[1]*KSSIG[2] - N[2]*KSSIG[1]; - NCSIG[1] = N[2]*KSSIG[0] - N[0]*KSSIG[2]; - NCSIG[2] = N[0]*KSSIG[1] - N[1]*KSSIG[0]; - - //VCS crossproduct of v[k] and KSS = v[k]XKSS - VCS[0] = v[1]*KSS[2] - v[2]*KSS[1]; - VCS[1] = v[2]*KSS[0] - v[0]*KSS[2]; - VCS[2] = v[0]*KSS[1] - v[1]*KSS[0]; - - //VCSIG crossproduct of v[k] and KSSIG = v[k]XKSSIG - VCSIG[0] = v[1]*KSSIG[2] - v[2]*KSSIG[1]; - VCSIG[1] = v[2]*KSSIG[0] - v[0]*KSSIG[2]; - VCSIG[2] = v[0]*KSSIG[1] - v[1]*KSSIG[0]; - - SDNCV = KSS[0]*NCV[0]+KSS[1]*NCV[1]+KSS[2]*NCV[2]; - - SIGDNCV = KSSIG[0]*NCV[0]+KSSIG[1]*NCV[1]+KSSIG[2]*NCV[2]; - - NDV = N[0]*v[0] + N[1]*v[1] + N[2]*v[2]; - - XS2 = XS[0]*XS[0]+XS[1]*XS[1]+XS[2]*XS[2]; - XA2 = XA[0]*XA[0]+XA[1]*XA[1]+XA[2]*XA[2]; - - NXA = N[0]*XA[0]+N[1]*XA[1]+N[2]*XA[2]; - NXS = N[0]*XS[0]+N[1]*XS[1]+N[2]*XS[2]; - - VDS = v[0]*KSS[0]+v[1]*KSS[1]+v[2]*KSS[2]; - VDSIG = v[0]*KSSIG[0]+v[1]*KSSIG[1]+v[2]*KSSIG[2]; - - NDS = N[0]*KSS[0]+N[1]*KSS[1]+N[2]*KSS[2]; - NDSIG = N[0]*KSSIG[0]+N[1]*KSSIG[1]+N[2]*KSSIG[2]; - - for(k=0;k<3;k++) - { - C1_5[k] = (N[k]*(12.0*SDNCV+6.0*DM*SIGDNCV/M)+9.0*NDV*NCS[k]+3.0*DM*NDV*NCSIG[k]/M -7.0*VCS[k]-3.0*DM*VCSIG[k]/M)/r3; - C2[k] = -MOR*MOR*MOR/r*3.0*eta*(N[k]*(XS2-XA2-5.0*NXS*NXS+5.0*NXA*NXA)+2.0*(XS[k]*NXS-XA[k]*NXA)); - C2_5[k] = (N[k]*(SDNCV*(-30.0*eta*NDV*NDV+24.0*eta*V1_V22-MOR*(38.0+25.0*eta))+DM/M*SIGDNCV*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22 - -MOR*(18.0+14.5*eta)))+NDV*v[k]*(SDNCV*(-9.0+9.0*eta)+DM/M* SIGDNCV*(-3.0+6.0*eta))+NCV[k]*(NDV*VDS*(-3.0+3.0*eta) - -8.0*MOR*eta*NDS-DM/M*(4.0*MOR*eta*NDSIG+3.0*NDV*VDSIG))+NDV*NCS[k]*(-22.5*eta*NDV*NDV+21.0*eta*V1_V22-MOR*(25.0+15.0*eta)) - +DM/M*NDV*NCSIG[k]*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22-MOR*(9.0+8.5*eta))+VCS[k]*(16.5*eta*NDV*NDV+MOR*(21.0+9.0*eta) - -14.0*eta*V1_V22)+DM/M*VCSIG[k]*(9.0*eta*NDV*NDV-7.0*eta*V1_V22+MOR*(9.0+4.5*eta)))/r3; - - if(Van_QM==1) - { - - if(m1>m2) - { - QMAux2_1[k] = (1.0-5.0*rS1*rS1)*N[k]+2.0*rS1*S1Dir[k]; - QMAux2_2[k] = (1.0-5.0*rS2*rS2)*N[k]+2.0*rS2*S2Dir[k]; - QMAux1[k] = Spin1Abs2*QMAux2_1[k]/nu+Spin2Abs2*QMAux2_2[k]*nu; - QM[k] = -1.5*MOR*MOR*MOR*eta*QMAux1[k]/r; - } - else - { - QMAux2_1[k] = (1.0-5.0*rS2*rS2)*N[k]+2.0*rS2*S2Dir[k]; - QMAux2_2[k] = (1.0-5.0*rS1*rS1)*N[k]+2.0*rS1*S1Dir[k]; - QMAux1[k] = Spin2Abs2*QMAux2_1[k]/nu+Spin1Abs2*QMAux2_2[k]*nu; - QM[k] = -1.5*MOR*MOR*MOR*eta*QMAux1[k]/r; - } - - } /* if(Van_QM==1) */ - - } /* k */ - - } /* if(Van_Spin==1) */ - - - - - -for(k=0;k<3;k++) - { - -if(usedOrNot[0] == 1) // PN0 (Newton) ~1/c^0 - { - a_pn1[0][k] = -m2*x[k]/r3; - a_pn2[0][k] = m1*x[k]/r3; - } - -if(usedOrNot[1] == 1) // PN1 ~1/c^2 - { - a_pn1[1][k] = ((AK2*N[k] + BK2*v[k])/r2)*m2; - a_pn2[1][k] = -((AK2*N[k] + BK2*v[k])/r2)*m1; - } - -if(usedOrNot[2] == 1) // PN2 ~1/c^4 - { - a_pn1[2][k] = ((AK4*N[k] + BK4*v[k])/r2)*m2; - a_pn2[2][k] = -((AK4*N[k] + BK4*v[k])/r2)*m1; - } - -if(usedOrNot[3] == 1) // PN2.5 ~1/c^5 - { - a_pn1[3][k] = ((AK5*N[k] + BK5*v[k])/r2)*m2; - a_pn2[3][k] = -((AK5*N[k] + BK5*v[k])/r2)*m1; - } - -if(usedOrNot[4] == 1) // PN3 ~1/c^6 - { - a_pn1[4][k] = ((AK6*N[k] + BK6*v[k])/r2)*m2; - a_pn2[4][k] = -((AK6*N[k] + BK6*v[k])/r2)*m1; - } - -if(usedOrNot[5] == 1) // PN3.5 ~1/c^7 - { - a_pn1[5][k] = ((AK7*N[k] + BK7*v[k])/r2)*m2; - a_pn2[5][k] = -((AK7*N[k] + BK7*v[k])/r2)*m1; - } - -if(Van_Spin == 1) // All the SPIN terms - { - a_pn1[6][k] += (C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4 + QM[k]/c_4)*m2/M; - a_pn2[6][k] += -(C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4 + QM[k]/c_4)*m1/M; - } - - A[k] = MOR*((AK2+AK4+AK5+AK6+AK7)*N[k] + (BK2+BK4+BK5+BK6+BK7)*v[k])/r + C1_5[k]/c_2 + C2[k]/c_4 + C2_5[k]/c_4 + QM[k]/c_4; - } - -// PN accelerations - - - - -// PN jerks - -for(k=0;k<3;k++) - { - AT[k] = A[k] - MOR*N[k]/r; // miert van AT - ? - } - -/* -AT[0] = A[0]; -AT[1] = A[1]; -AT[2] = A[2]; -*/ - -RPP = V1_V22/r + AT[0]*N[0]+AT[1]*N[1] + AT[2]*N[2] - RP*RP/r; -VA = AT[0]*v[0] + AT[1]*v[1] + AT[2]*v[2]; - -for(k=0;k<3;k++) NDOT[k] = (v[k]-N[k]*RP)/r; - -NVDOT = NDOT[0]*v[0]+NDOT[1]*v[1]+NDOT[2]*v[2]+N[0]*AT[0]+N[1]*AT[1]+N[2]*AT[2]; - -//NDOTCV crossproduct of NDOT[k] and relative v = NDOT[k]Xv[j] -NDOTCV[0] = NDOT[1]*v[2] - NDOT[2]*v[1]; -NDOTCV[1] = NDOT[2]*v[0] - NDOT[0]*v[2]; -NDOTCV[2] = NDOT[0]*v[1] - NDOT[1]*v[0]; - -//NCA crossproduct of N and AT = N[k]XAT[j] -NCA[0] = N[1]*AT[2] - N[2]*AT[1]; -NCA[1] = N[2]*AT[0] - N[0]*AT[2]; -NCA[2] = N[0]*AT[1] - N[1]*AT[0]; - -ADK2 = 0.0; BDK2 = 0.0; -ADK4 = 0.0; BDK4 = 0.0; -ADK5 = 0.0; BDK5 = 0.0; -ADK6 = 0.0; BDK6 = 0.0; -ADK7 = 0.0; BDK7 = 0.0; - - -for(k=0;k<3;k++) - { - C1_5D[k] = 0.0; - C2D[k] = 0.0; - C2_5D[k] = 0.0; - QMD[k] = 0.0; - } - -if(usedOrNot[1] == 1) // PN1 ~1/c^2 - { - A1D = -2.0*(2.0+eta)*MOR*RP/r - 2.0*(1.0+3.0*eta)*VA + 3.0*eta*RP*RPP; - B1D = 2.0*(2.0-eta)*RPP; - - ADK2 = A1D/c_2; - BDK2 = B1D/c_2; - } - -if(usedOrNot[2] == 1) // PN2 ~1/c^4 - { - A2D = 1.5*(12.0+29.0*eta)*MOR*MOR*RP/r -eta*(3.0-4.0*eta)*4.0*V1_V22*VA - 7.5*eta*(1.0-3.0*eta)*RPP -0.5*eta*(13.0-4.0*eta)*MOR*RP*V1_V22/r+eta*(13.0-4.0*eta)*MOR*VA -(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RP*RP/r+2.0*(2.0+25.0*eta+2.0*eta*eta)*MOR*RP*RPP + 3.0*eta*(3.0-4.0*eta)*VA*RP*RP + 3.0*eta*(3.0-4.0*eta)*V1_V22*RP*RPP; - B2D = -0.5*RPP*((4.0+41.0*eta+8.0*eta*eta)*MOR - eta*(15.0+4.0*eta)*V1_V22+3.0*eta*(3.0+2.0*eta)*RP*RP) - 0.5*RP*(-(4.0+41.0*eta+8.0*eta*eta)*MOR*RP/r - 2.0*eta*(15.0+4.0*eta)*VA + 6.0*eta*(3.0+2.0*eta)*RP*RPP); - - ADK4 = A2D/c_4; - BDK4 = B2D/c_4; - } - -if(usedOrNot[3] == 1) // PN2.5 ~1/c^5 - { - A2_5D = -1.6*eta*MOR*RP*RP*(17.0/3.0*MOR+3.0*V1_V22)/r +1.6*eta*MOR*RPP*(17.0/3.0*MOR+3.0*V1_V22)+1.6*eta*MOR*RP*(-17.0*MOR*RP/3.0/r+6.0*VA); - B2_5D = 1.6*eta*MOR*RP*(3.0*MOR+V1_V22)/r - 1.6*eta*MOR*(-3.0*MOR*RP/r+2.0*VA); - - ADK5 = A2_5D/c_5; - BDK5 = B2_5D/c_5; - } - -if(usedOrNot[4] == 1) // PN3 ~1/c^6 - { - A3D = 6.0*eta*RP*RP*RP*RP*RP*RPP*(35.0-175.0*eta+175.0*eta*eta)/16.0 + eta*(4.0*RP*RP*RP*RPP*V1_V22 + 2.0*RP*RP*RP*RP*VA)*(-15.0+135.0*eta/2.0-255.0*eta*eta/4.0)/2.0 + eta*(2.0*RP*RPP*V1_V22*V1_V22+4.0*RP*RP*V1_V22*VA)/2.0*(15.0-237.0*eta/2.0+45.0*eta*eta) + 6.0*V1_V22*V1_V22*VA*eta*(-11.0/4.0-49.0*eta/4.0-13.0*eta*eta) + MOR*(4.0*RP*RP*RP*RPP*eta*(-79.0+69.0/2.0*eta+30.0*eta*eta) + eta*(2.0*RP*RPP*V1_V22+2.0*RP*RP*VA)*(121.0-16.0*eta-20.0*eta*eta)+4.0*V1_V22*VA*eta*(-75.0/4.0-8.0*eta+10.0*eta*eta)) - MOR*RP*((-79.0+69.0*eta/2.0+30.0*eta*eta)*RP*RP*RP*RP*eta+eta*RP*RP*V1_V22*(121.0-16.0*eta-20.0*eta*eta)+eta*V1_V22*V1_V22*(-75.0/4.0-8.0*eta+10.0*eta*eta))/r - 2.0*MOR*MOR*RP*(RP*RP*((-1.0-615.0*PI2*eta/64.0)-22717.0*eta/168.0-11.0*eta*eta/8.0+7.0*eta*eta*eta)+eta*V1_V22*((20827.0/840.0+123.0*PI2/64.0)-eta*eta))/r + MOR*MOR*(2.0*RP*RPP*((-1.0-615*PI2*eta/64.0)-22717.0*eta/168.0-11.0*eta*eta/8.0+7*eta*eta*eta)+2.0*eta*VA*((20827.0/840.0 +123.0*PI2/64.0)-eta*eta)) - 3.0*MOR*MOR*MOR*RP*(16.0+(1399.0/12.0-41.0*PI2/16.0)*eta+71.0*eta*eta/2.0)/r; - B3D = 75.0*RP*RP*RP*RP*RPP*eta*(3.0/8.0-eta-.25*eta*eta)+eta*(3.0*RP*RP*RPP*V1_V22+2.0*RP*RP*RP*VA)*(-12.0+111.0*eta/4.0+12.0*eta*eta)+eta*(RPP*V1_V22*V1_V22+4.0*RP*V1_V22*VA)*(65.0/8.0-19.0*eta-6.0*eta*eta)-MOR*RP*(RP*RP*RP*eta*(-329.0/6.0-59.0*eta/2.0-18.0*eta*eta)+RP*V1_V22*eta*(15.0+27.0*eta+10.0*eta*eta))/r+MOR*(3.0*RP*RP*RPP*eta*(-329.0/6.0-59.0*eta/2.0-18.0*eta*eta)+eta*(RPP*V1_V22+2.0*RP*VA)*(15.0+27.0*eta+10.0*eta*eta))-2.0*MOR*MOR*RP*(RP*((4.0+123.0*PI2*eta/32.0)+5849.0*eta/840.0-25.0*eta*eta-8.0*eta*eta*eta))/r+MOR*MOR*(RPP*((4.0+123.0*PI2*eta/32.0)+5849.0/840.0*eta-25.0*eta*eta-8.0*eta*eta*eta)); - - ADK6 = A3D/c_6; - BDK6 = B3D/c_6; - } - -if(usedOrNot[5] == 1) // PN3.5 ~1/c^7 - { - A3_5D = MOR*eta*(-RP*(V1_V22*V1_V22*(-366.0/35.0-12.0*eta)+V1_V22*RP*RP*(114.0+12.0*eta)+RP*RP*RP*RP*(-112.0))/r+4.0*V1_V22*VA*(-366.0/35.0-12.0*eta)+2.0*(VA*RP*RP+RP*RPP*V1_V22)*(114.0+12.0*eta)+4.0*RP*RP*RP*RPP*(-112.0)+MOR*(2.0*VA*(-692.0/35.0+724.0*eta/15.0)+2.0*RP*RPP*(-294.0/5.0-376.0*eta/5.0)-2.0*RP*(V1_V22*(-692.0/35.0+724.0*eta/15.0)+RP*RP*(-294.0/5.0-376.0*eta/5.0))/r-3.0*MOR*RP*(-3956.0/35.0-184.0*eta/5.0)/r)); - B3_5D = MOR*eta*(4.0*V1_V22*VA*(626.0/35.0+12.0*eta/5.0)+2.0*(VA*RP*RP+V1_V22*RP*RPP)*(-678.0/5.0-12.0*eta/5.0)+4.0*RP*RP*RP*RPP*120.0-RP*(V1_V22*V1_V22*(626.0/35.0+12.0*eta/5.0)+V1_V22*RP*RP*(-678.0/5.0-12.0*eta/5.0)+120.0*RP*RP*RP*RP)/r+MOR*(2.0*VA*(-164.0/21.0-148.0*eta/5.0)+2*RP*RPP*(82.0/3.0+848.0*eta/15.0)-2.0*RP*(V1_V22*(-164.0/21-148.0*eta/5.0)+RP*RP*(82.0/3.0+848.0*eta/15.0))/r-3.0*MOR*RP*(1060.0/21.0+104.0*eta/5.0)/r)); - - ADK7 = A3_5D/c_7; - BDK7 = B3_5D/c_7; - } - - - - - - if(Van_Spin==1) - { - - //L crossproduct of x[k] and relative v = x[k]Xv[j] - L[0] = x[1]*v[2] - x[2]*v[1]; - L[1] = x[2]*v[0] - x[0]*v[2]; - L[2] = x[0]*v[1] - x[1]*v[0]; - - LABS = sqrt(L[0]*L[0]+L[1]*L[1]+L[2]*L[2]); - - LU[0] = L[0]/LABS; - LU[1] = L[1]/LABS; - LU[2] = L[2]/LABS; - - S1DLU = S1[0]*LU[0]+S1[1]*LU[1]+S1[2]*LU[2]; - S2DLU = S2[0]*LU[0]+S2[1]*LU[1]+S2[2]*LU[2]; - - for(k=0;k<3;k++) - { - SU1[k] = MOR*eta*(N[k]*(-4.0*VDS-2.0*DM/M*VDSIG)+ v[k]*(3.0*NDS+DM/M*NDSIG)+NDV*(2.0*KSS[k]+DM/M*KSSIG[k])) /r; - SV1[k] = MOR*(N[k]*(VDSIG*(-2.0+4.0*eta)-2.0*DM/M*VDS)+ v[k]*(NDSIG*(1.0-eta)+DM/M*NDS)+NDV*(KSSIG[k]*(1.0- 2.0*eta)+ DM/M*KSS[k]))/r; - - SS1[k] = 0.5*(L[k]*(4.0+3.0*(m2/m1))+ (S2[k]-3.0*S2DLU*LU[k]))/r3; - SS2[k] = 0.5*(L[k]*(4.0+3.0*(m1/m2))+ (S1[k]-3.0*S1DLU*LU[k]))/r3; - - SU2[k] = MOR*eta/r*(N[k]*(VDS*(-2.0*V1_V22+3.0*NDV*NDV- 6.0*eta*NDV*NDV+7.0*MOR-8.0*eta*MOR)-14.0*MOR*NDS*NDV+ DM/M*VDSIG*eta*(-3.0*NDV*NDV-4.0*MOR)+DM/M*MOR*NDSIG*NDV* (2.0-eta/2.))+v[k]*(NDS*(2.0*V1_V22-4.0*eta*V1_V22-3.0*NDV* NDV+7.5*eta*NDV*NDV+4.0*MOR-6.0*eta*MOR)+VDS*NDV*(2.0- 6.0*eta)+ DM/M*NDSIG*(-1.5*eta*V1_V22+3.0*eta*NDV*NDV-MOR-3.5*eta* MOR)-3.0*DM/M*VDSIG*NDV*eta)+KSS[k]*NDV*(V1_V22-2.0*eta* V1_V22-1.5*NDV*NDV+3.0*eta*NDV*NDV-MOR+2.0*eta*MOR)+ DM/M*KSSIG[k]*NDV*(-eta*V1_V22+1.5*eta*NDV*NDV+ (eta-1.)*MOR)); - SV2[k] = MOR/r*(N[k]*(VDSIG*eta*(-2.0*V1_V22+6.0*eta*NDV* NDV+(3.0+8.0*eta)*MOR)+MOR*NDSIG*NDV*(2.0-22.5*eta+2.0* eta*eta)+ DM/M*VDS*eta*(-3.0*NDV*NDV-4.0*MOR)+DM/M*MOR*NDS*NDV*(2.0- 0.5*eta))+v[k]*(NDSIG*(0.5*eta*V1_V22+2.0*eta*eta*V1_V22- 4.5*eta*eta*NDV*NDV+(4.5*eta-1.0+8.0*eta*eta)*MOR)+VDSIG*NDV* eta*(6.0*eta-1.)-3.0*DM/M*VDS*NDV*eta+DM/M*NDS*(-1.5* eta*V1_V22+ 3.0*eta*NDV*NDV-(1.0+3.5*eta)*MOR))+KSSIG[k]*NDV*(2.0*eta*eta* V1_V22-3.0*eta*eta*NDV*NDV+(-1.0+4.0*eta-2.0*eta*eta)*MOR)+ DM/M*KSS[k]*NDV*(-eta*V1_V22+1.5*eta*NDV*NDV+(-1.0+eta)* MOR)); - } - - //SS1 crossproduct of SS1 and S1 = SS1[k]XS1[j] - SS1aux[0] = SS1[1]*S1[2] - SS1[2]*S1[1]; - SS1aux[1] = SS1[2]*S1[0] - SS1[0]*S1[2]; - SS1aux[2] = SS1[0]*S1[1] - SS1[1]*S1[0]; - - SS1[0] = SS1aux[0]; - SS1[1] = SS1aux[1]; - SS1[2] = SS1aux[2]; - - //SS2 crossproduct of SS2 and S2 = SS2[k]XS2[j] - SS2aux[0] = SS2[1]*S2[2] - SS2[2]*S2[1]; - SS2aux[1] = SS2[2]*S2[0] - SS2[0]*S2[2]; - SS2aux[2] = SS2[0]*S2[1] - SS2[1]*S2[0]; - - SS2[0] = SS2aux[0]; - SS2[1] = SS2aux[1]; - SS2[2] = SS2aux[2]; - - SPINPrev[0][0] = SPIN[0][0]; - SPINPrev[1][0] = SPIN[1][0]; - SPINPrev[2][0] = SPIN[2][0]; - - SPINPrev[0][1] = SPIN[0][1]; - SPINPrev[1][1] = SPIN[1][1]; - SPINPrev[2][1] = SPIN[2][1]; - - SpinPrev2_1 = SPINPrev[0][0]*SPINPrev[0][0] + SPINPrev[1][0]*SPINPrev[1][0] + SPINPrev[2][0]*SPINPrev[2][0]; - SpinPrev2_2 = SPINPrev[0][1]*SPINPrev[0][1] + SPINPrev[1][1]*SPINPrev[1][1] + SPINPrev[2][1]*SPINPrev[2][1]; - - SPSPP1 = 0.0; - SPSPP2 = 0.0; - - Spin1AbsNew2 = 0.0; - Spin2AbsNew2 = 0.0; - - for(k=0;k<3;k++) - { - SU[k] = SU1[k]/c_2 + SU2[k]/c_4 + (SS1[k] + SS2[k])/c_2; - SV[k] = SV1[k]/c_2 + SV2[k]/c_4+M*(SS2[k]/m2-SS1[k]/ m1)/c_2; - - KSS[k] = KSS[k] + SU[k]*dt_bh; // integrate for dt_bh timestep - KSSIG[k] = KSSIG[k] + SV[k]*dt_bh; - - SPIN[k][0] = m1*(M*KSS[k]-m2*KSSIG[k])/M/M/m1/m1*c_1; - SPIN[k][1] = m2*(M*KSS[k]+m1*KSSIG[k])/M/M/m2/m2*c_1; - Spin1AbsNew2 += SPIN[k][0]*SPIN[k][0]; - Spin2AbsNew2 += SPIN[k][1]*SPIN[k][1]; - XAD[k] = 0.5/(M*M*m1*m2)*(-SU[k]*M*DM-SV[k]*(m1*m1+m2*m2)); - XSD[k] = 0.5/(M*M*m1*m2)*(SU[k]*M*M+SV[k]*(m1*m1-m2*m2)); - - if(m1>m2) - { - SPSPP1 += SPINPrev[k][0]*(SPIN[k][0]-SPINPrev[k][0])/dt_bh; - SPSPP2 += SPINPrev[k][1]*(SPIN[k][1]-SPINPrev[k][1])/dt_bh; - } - else - { - SPSPP1 += SPINPrev[k][1]*(SPIN[k][1]-SPINPrev[k][1])/dt_bh; - SPSPP2 += SPINPrev[k][0]*(SPIN[k][0]-SPINPrev[k][0])/dt_bh; - } - } - - Spin1AbsNew = sqrt(Spin1AbsNew2); - Spin2AbsNew = sqrt(Spin2AbsNew2); - - for(k=0;k<3;k++) - { - S1DirNew[k] = SPIN[k][0]/Spin1AbsNew; - S2DirNew[k] = SPIN[k][1]/Spin2AbsNew; - } - - - //NDOTCS crossproduct of NDOT and KSS = NDOT[k]XKSS[j] - NDOTCS[0] = NDOT[1]*KSS[2] - NDOT[2]*KSS[1]; - NDOTCS[1] = NDOT[2]*KSS[0] - NDOT[0]*KSS[2]; - NDOTCS[2] = NDOT[0]*KSS[1] - NDOT[1]*KSS[0]; - //NCSU crossproduct of N and SU = N[k]XSU[j] - NCSU[0] = N[1]*SU[2] - N[2]*SU[1]; - NCSU[1] = N[2]*SU[0] - N[0]*SU[2]; - NCSU[2] = N[0]*SU[1] - N[1]*SU[0]; - //NDOTCSIG crossproduct of NDOT and KSSIG = NDOT[k]XKSSIG[j] - NDOTCSIG[0] = NDOT[1]*KSSIG[2] - NDOT[2]*KSSIG[1]; - NDOTCSIG[1] = NDOT[2]*KSSIG[0] - NDOT[0]*KSSIG[2]; - NDOTCSIG[2] = NDOT[0]*KSSIG[1] - NDOT[1]*KSSIG[0]; - //NCSV crossproduct of N and SV = N[k]XSV[j] - NCSV[0] = N[1]*SV[2] - N[2]*SV[1]; - NCSV[1] = N[2]*SV[0] - N[0]*SV[2]; - NCSV[2] = N[0]*SV[1] - N[1]*SV[0]; - //ACS crossproduct of AT and KSS = AT[k]XKSS[j] - ACS[0] = AT[1]*KSS[2] - AT[2]*KSS[1]; - ACS[1] = AT[2]*KSS[0] - AT[0]*KSS[2]; - ACS[2] = AT[0]*KSS[1] - AT[1]*KSS[0]; - //VCSU crossproduct of relative v and SU = v[k]XSU[j] - VCSU[0] = v[1]*SU[2] - v[2]*SU[1]; - VCSU[1] = v[2]*SU[0] - v[0]*SU[2]; - VCSU[2] = v[0]*SU[1] - v[1]*SU[0]; - //ACSIG crossproduct of AT and KSSIG = AT[k]XKSSIG[j] - ACSIG[0] = AT[1]*KSSIG[2] - AT[2]*KSSIG[1]; - ACSIG[1] = AT[2]*KSSIG[0] - AT[0]*KSSIG[2]; - ACSIG[2] = AT[0]*KSSIG[1] - AT[1]*KSSIG[0]; - //VCSV crossproduct of relative v and SV = v[k]XSV[j] - VCSV[0] = v[1]*SV[2] - v[2]*SV[1]; - VCSV[1] = v[2]*SV[0] - v[0]*SV[2]; - VCSV[2] = v[0]*SV[1] - v[1]*SV[0]; - - SNVDOT = SU[0]*NCV[0]+SU[1]*NCV[1]+SU[2]*NCV[2]+ KSS[0]*NDOTCV[0]+KSS[1]*NDOTCV[1]+KSS[2]*NDOTCV[2]+ KSS[0]*NCA[0]+KSS[1]*NCA[1]+KSS[2]*NCA[2]; - - SIGNVDOT = SV[0]*NCV[0]+SV[1]*NCV[1]+SV[2]*NCV[2]+ KSSIG[0]*NDOTCV[0]+KSSIG[1]*NDOTCV[1]+KSSIG[2]*NDOTCV[2]+ KSSIG[0]*NCA[0]+KSSIG[1]*NCA[1]+KSSIG[2]*NCA[2]; - - NSDOT = NDOT[0]*KSS[0]+NDOT[1]*KSS[1]+NDOT[2]*KSS[2]+ N[0]*SU[0]+N[1]*SU[1]+N[2]*SU[2]; - NSIGDOT = NDOT[0]*KSSIG[0]+NDOT[1]*KSSIG[1]+NDOT[2]*KSSIG[2]+ N[0]*SV[0]+N[1]*SV[1]+N[2]*SV[2]; - VSDOT = AT[0]*KSS[0]+AT[1]*KSS[1]+AT[2]*KSS[2]+ v[0]*SU[0]+v[1]*SU[1]+v[2]*SU[2]; - VSIGDOT = AT[0]*KSSIG[0]+AT[1]*KSSIG[1]+AT[2]*KSSIG[2]+ v[0]*SV[0]+v[1]*SV[1]+v[2]*SV[2]; - - NXSDOT = NDOT[0]*XS[0]+NDOT[1]*XS[1]+NDOT[2]*XS[2]+ N[0]*XSD[0]+N[1]*XSD[1]+N[2]*XSD[2]; - NXADOT = NDOT[0]*XA[0]+NDOT[1]*XA[1]+NDOT[2]*XA[2]+ N[0]*XAD[0]+N[1]*XAD[1]+N[2]*XAD[2]; - - rS1p = -rS1*NDV/r; - rS2p = -rS2*NDV/r; - - for(k=0;k<3;k++) - { - S1p[k] = (S1DirNew[k] - S1Dir[k])/dt_bh; - S2p[k] = (S2DirNew[k] - S2Dir[k])/dt_bh; - - rS1p += v[k]*S1Dir[k]/r + N[k]*S1p[k]; - rS2p += v[k]*S2Dir[k]/r + N[k]*S2p[k]; - - Np[k] = (v[k] - N[k]*NDV)/r; - } - - for(k=0;k<3;k++) - { - C1_5D[k] = -3.0*RP/r*C1_5[k]+(NDOT[k]*(12.0*SDNCV+6.0*DM/M* SIGDNCV)+N[k]*(12.0*SNVDOT+6.0*DM/M*SIGNVDOT)+9.0*NVDOT* NCS[k]+9.0*NDV*(NDOTCS[k]+NCSU[k])+3.0*DM/M*(NVDOT*NCSIG[k]+ NDV*(NDOTCSIG[k]+NCSV[k]))-7.0*(ACS[k]+VCSU[k])-3.0*DM/M* (ACSIG[k]+VCSV[k]))/(r3); - C2D[k] = -4.0*RP/r*C2[k]-MOR*MOR*MOR*3.0*eta/r*(NDOT[k]* (XS2-XA2-5.0*NXS*NXS+5.0*NXA*NXA)+N[k]*(2.0*(XS[0]*XSD[0]+ XS[1]*XSD[1]+XS[2]*XSD[2]-XA[0]*XAD[0]-XA[1]*XAD[1]- XA[2]*XAD[2])-10.0*NXS*NXSDOT+10.0*NXA*NXADOT)+2.0*(XSD[k]* NXS+XS[k]*NXSDOT-XAD[k]*NXA-XA[k]*NXADOT)); - C2_5D[k] = -3.0*RP/r*C2_5[k]+(NDOT[k]*(SDNCV*(-30.0*eta* NDV*NDV+24.0*eta*V1_V22-MOR*(38.0+25.0*eta))+DM/M*SIGDNCV* (-15.0*eta*NDV*NDV+12.0*eta*V1_V22-MOR*(18.0+14.5*eta)))+ N[k]*(SNVDOT*(-30.0*eta*NDV*NDV+24.0*eta*V1_V22-MOR* (38.0+25.0*eta))+SDNCV*(-60.0*eta*NDV*NVDOT+48.0*eta*VA+ MOR*RP/r*(38.0+25.0*eta))+DM/M*SIGNVDOT*(-15.0*eta*NDV* NDV+12.0*eta*V1_V22-MOR*(18.0+14.5*eta))+DM/M*SIGDNCV* (-30.0*eta*NDV*NVDOT+24.0*eta*VA+MOR*RP/r*(18.0+14.5*eta)))+ (NVDOT*v[k]+NDV*AT[k])*(SDNCV*(-9.0+9.0*eta)+DM/M*SIGDNCV* (-3.0+6.0*eta))+NDV*v[k]*(SNVDOT*(-9.0+9.0*eta)+DM/M* SIGNVDOT*(-3.0+6.0*eta))+(NDOTCV[k]+NCA[k])*(NDV*VDS*(-3.0+ 3.0*eta)-8.0*MOR*eta*NDS-DM/M*(4.0*MOR*eta*NDSIG+3.0*NDV*VDSIG) )+NCV[k]*((NVDOT*VDS+NDV*VSDOT)*(-3.0+3.0*eta)-8.0*eta*MOR* (NSDOT-RP/r*NDS)-DM/M*(4.0*eta*MOR*(NSIGDOT-RP/r*NDSIG)+ 3.0*(NVDOT*VDSIG+NDV*VSIGDOT)))+(NVDOT*NCS[k]+NDV* (NDOTCS[k]+NCSU[k]))*(-22.5*eta*NDV*NDV+21.0*eta*V1_V22- MOR*(25.0+15.0*eta))+NDV*NCS[k]*(-45.0*eta*NDV*NVDOT+42.0*eta* VA+MOR*RP/r*(25.0+15.0*eta))+DM/M*(NVDOT*NCSIG[k]+NDV* (NDOTCSIG[k]+NCSV[k]))*(-15.0*eta*NDV*NDV+12.0*eta*V1_V22- MOR*(9.0+8.5*eta))+DM/M*NDV*NCSIG[k]*(-30.0*eta*NDV*NVDOT+ 24.0*eta*VA+MOR*RP/r*(9.0+8.5*eta))+(ACS[k]+VCSU[k])* (16.5*eta*NDV*NDV+MOR*(21.0+9.0*eta)-14.0*eta*V1_V22)+ VCS[k]*(33.0*eta*NDV*NVDOT-MOR*RP/r*(21.0+9.0*eta)- 28.0*eta*VA)+DM/M*(ACSIG[k]+VCSV[k])*(9.0*eta*NDV*NDV- 7.0*eta*V1_V22+MOR*(9.0+4.5*eta))+DM/M*VCSIG[k]*(18.0* eta*NDV*NVDOT-14.0*eta*VA-MOR*RP/r*(9.0+4.5*eta)))/ (r3); - - - if(Van_QM==1) - { - - if(m1>m2) - { - QMD[k] = -1.5*MOR*MOR*MOR*eta*(-4.0*RP*QMAux1[k]/r2+( 2.0*(SPSPP1*QMAux2_1[k]/nu+SPSPP2*QMAux2_2[k]*nu) + SpinPrev2_1*(-10.0*rS1*rS1p*N[k]+(1.0-5.0*rS1*rS1)*Np[k]+2.0*rS1p*S1Dir[k]+2.0*rS1*S1p[k])/nu + SpinPrev2_2*(-10.0*rS2*rS2p*N[k]+(1.0-5.0*rS2*rS2)*Np[k]+2.0*rS2p*S2Dir[k]+2.0*rS2*S2p[k])*nu )/r); - } - else - { - QMD[k] = -1.5*MOR*MOR*MOR*eta*(-4.0*RP*QMAux1[k]/r2+( 2.0*(SPSPP2*QMAux2_1[k]/nu+SPSPP1*QMAux2_2[k]*nu) + SpinPrev2_2*(-10.0*rS2*rS2p*N[k]+(1.0-5.0*rS2*rS2)*Np[k]+2.0*rS2p*S2Dir[k]+2.0*rS2*S2p[k])/nu + SpinPrev2_1*(-10.0*rS1*rS1p*N[k]+(1.0-5.0*rS1*rS1)*Np[k]+2.0*rS1p*S1Dir[k]+2.0*rS1*S1p[k])*nu )/r); - } - - } /* if(Van_QM==1) */ - - } /* k */ - - - } /* if(Van_Spin==1) */ - - - ADK = ADK2+ADK4+ADK5+ADK6+ADK7; - BDK = BDK2+BDK4+BDK5+BDK6+BDK7; - - KSAK = AK2+AK4+AK5+AK6+AK7; - KSBK = BK2+BK4+BK5+BK6+BK7; - - for(k=0;k<3;k++) AD[k] = -2.0*MOR*RP*(KSAK*N[k]+KSBK*v[k])/r2 + MOR*(ADK*N[k]+BDK*v[k])/r + MOR*(KSAK*(v[k]-N[k]*RP)/r+KSBK*AT[k])/r + C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4 + QMD[k]/c_4; - - -for(k=0;k<3;k++) // new values of the BH's spins, returned back to the main program... - { - spin1[k] = SPIN[k][0]; - spin2[k] = SPIN[k][1]; - } - - - - - - - -for(k=0;k<3;k++) - { - -if(usedOrNot[0] == 1) // PN0 (Newton) ~1/c^0 - { - adot_pn1[0][k] = -m2*(v[k]/r3 - 3.0*RP*x[k]/r2/r2); - adot_pn2[0][k] = m1*(v[k]/r3 - 3.0*RP*x[k]/r2/r2); - } - -if(usedOrNot[1] == 1) // PN1 ~1/c^2 - { - adot_pn1[1][k] = (-2.0*MOR*RP*(AK2*N[k]+BK2*v[k])/r2 + MOR*(ADK2*N[k]+BDK2*v[k])/r + MOR*(AK2*(v[k]-N[k]*RP)/r+BK2*A[k])/r)*m2/M; - adot_pn2[1][k] = -(-2.0*MOR*RP*(AK2*N[k]+BK2*v[k])/r2 + MOR*(ADK2*N[k]+BDK2*v[k])/r + MOR*(AK2*(v[k]-N[k]*RP)/r+BK2*A[k])/r)*m1/M; - } - -if(usedOrNot[2] == 1) // PN2 ~1/c^4 - { - adot_pn1[2][k] = (-2.0*MOR*RP*(AK4*N[k]+BK4*v[k])/r2 + MOR*(ADK4*N[k]+BDK4*v[k])/r + MOR*(AK4*(v[k]-N[k]*RP)/r+BK4*A[k])/r)*m2/M; - adot_pn2[2][k] = -(-2.0*MOR*RP*(AK4*N[k]+BK4*v[k])/r2 + MOR*(ADK4*N[k]+BDK4*v[k])/r + MOR*(AK4*(v[k]-N[k]*RP)/r+BK4*A[k])/r)*m1/M; - } - -if(usedOrNot[3] == 1) // PN2.5 ~1/c^5 - { - adot_pn1[3][k] = (-2.0*MOR*RP*(AK5*N[k]+BK5*v[k])/r2 + MOR*(ADK5*N[k]+BDK5*v[k])/r + MOR*(AK5*(v[k]-N[k]*RP)/r+BK5*A[k])/r)*m2/M; - adot_pn2[3][k] = -(-2.0*MOR*RP*(AK5*N[k]+BK5*v[k])/r2 + MOR*(ADK5*N[k]+BDK5*v[k])/r + MOR*(AK5*(v[k]-N[k]*RP)/r+BK5*A[k])/r)*m1/M; - } - -if(usedOrNot[4] == 1) // PN3 ~1/c^6 - { - adot_pn1[4][k] = (-2.0*MOR*RP*(AK6*N[k]+BK6*v[k])/r2 + MOR*(ADK6*N[k]+BDK6*v[k])/r + MOR*(AK6*(v[k]-N[k]*RP)/r+BK6*A[k])/r)*m2/M; - adot_pn2[4][k] = -(-2.0*MOR*RP*(AK6*N[k]+BK6*v[k])/r2 + MOR*(ADK6*N[k]+BDK6*v[k])/r + MOR*(AK6*(v[k]-N[k]*RP)/r+BK6*A[k])/r)*m1/M; - } - -if(usedOrNot[5] == 1) // PN3.5 ~1/c^7 - { - adot_pn1[5][k] = (-2.0*MOR*RP*(AK7*N[k]+BK7*v[k])/r2 + MOR*(ADK7*N[k]+BDK7*v[k])/r + MOR*(AK7*(v[k]-N[k]*RP)/r+BK7*A[k])/r)*m2/M; - adot_pn2[5][k] = -(-2.0*MOR*RP*(AK7*N[k]+BK7*v[k])/r2 + MOR*(ADK7*N[k]+BDK7*v[k])/r + MOR*(AK7*(v[k]-N[k]*RP)/r+BK7*A[k])/r)*m1/M; - } - - -if(Van_Spin == 1) // All the SPIN terms - { - adot_pn1[6][k] += (C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4 + QMD[k]/c_4)*m2/M; - adot_pn2[6][k] += -(C1_5D[k]/c_2 + C2D[k]/c_4 +C2_5D[k]/c_4 + QMD[k]/c_4)*m1/M; - } - - - } - -// PN jerks - - - - -// Check RS_DIST conditions !!! - -RS_DIST = 4.0*(2.0*m1/c_2 + 2.0*m2/c_2); - -if(r < RS_DIST) - { - if(myRank == rootRank) - { - fprintf(stdout,"PN RSDIST: r = %.8E \t RS = %.8E \n", r, RS_DIST); - fflush(stdout); - } - return(505); - } -else - { - return(0); - } - - -} -/***************************************************************************/ diff --git a/star_destr.c b/star_destr.c deleted file mode 100644 index 612d765..0000000 --- a/star_destr.c +++ /dev/null @@ -1,275 +0,0 @@ -/***************************************************************************** - File Name : "Star Destr.c" - Contents : star "destruction" by tidal field of "live" BH (1 or 2) - Coded by : Peter Berczik - Last redaction : 2010.IX.14 1:43PM -*****************************************************************************/ - -void star_destr(double time, - int n, - int ind[], - double m[], - double x[][3], - double v[][3], - double pot[], - double a[][3], - double adot[][3], - double t[], - double dt[], - int N, - double m_N[], - double x_N[][3], - double v_N[][3], - double a_N[][3], - double adot_N[][3], - double t_N[], - double *m_bh, - int *num_bh, - int i_bh) -{ - -int n_end, k_act, N_end; - -double R_t, e_kin=0.0, e_pot_BH=0.0, e_pot=0.0, e_corr=0.0; - -double eps_bh, eps2, eps_bh2, rsb, rkb2, rks2, xp[3], v_bh[3]; - -//double vir = 0.6, gamma = 1000.0; - -double x_max = 1.0E+03, v_max = 1.0E-08, - a_max = 1.0E-08, adot_max = 1.0E-08; - - -if( time < t_diss_on ) return; - - -eps_bh = eps; - -eps2 = SQR(eps); -eps_bh2 = SQR(eps_bh); - -/* -m_s = 1.0/N; -R_s = 2.52E-08/2.507328103; -R_t = gamma * R_s * pow( 2.0*(*m_bh/m_s), over3 ); -*/ - -R_t = R_TIDAL; - -/* -if(myRank == rootRank) - { - printf("%.6E \t %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, i_bh, ind[i_bh], R_t, *m_bh, *num_bh); - fflush(stdout); - } -*/ - -#ifdef ADD_BH2 - n_end = n-2; - N_end = N-2; -#else -#ifdef ADD_BH1 - n_end = n-1; - N_end = N-1; -#endif // ADD_BH1 -#endif // ADD_BH2 - -/* -if(myRank == rootRank) - { - printf("%.6E \t %06d %06d %06d %06d \t %.6E \t %.6E %06d \n", time, n, n_end, i_bh, ind[i_bh], R_t, *m_bh, *num_bh); - fflush(stdout); - } -*/ - -for(i=0; i