first commit
This commit is contained in:
commit
aaf40315c8
15 changed files with 17524 additions and 0 deletions
30
Makefile
Normal file
30
Makefile
Normal file
|
|
@ -0,0 +1,30 @@
|
||||||
|
CUDAHOME ?= /usr/local/cuda
|
||||||
|
CPPFLAGS += -DYEBISU -DETICS
|
||||||
|
OPTIMIZATION ?= 3
|
||||||
|
|
||||||
|
DTSCF ?= 0.015625
|
||||||
|
|
||||||
|
CUDAINC = -I$(CUDAHOME)/include -I$(CUDAHOME)/samples/common/inc/
|
||||||
|
CUDALIB = -L$(CUDAHOME)/lib64 -lcudart -lcudadevrt
|
||||||
|
|
||||||
|
GRAPEHOME = ../grapite
|
||||||
|
GRAPELIB = -L$(GRAPEHOME) -lgrapite
|
||||||
|
yebisu: GRAPEHOME = ../yebisu
|
||||||
|
yebisu: GRAPELIB = -L$(GRAPEHOME) -lyebisug6
|
||||||
|
GRAPEINC = -I$(GRAPEHOME)
|
||||||
|
|
||||||
|
CFLAGS ?= -mcmodel=large
|
||||||
|
CFLAGS += -O$(OPTIMIZATION)
|
||||||
|
INC = $(GRAPEINC) $(CUDAINC)
|
||||||
|
LIB = $(GRAPELIB) $(CUDALIB) -lm -lgcc -lgfortran -lstdc++
|
||||||
|
MPICC ?= mpicc
|
||||||
|
EXECUTABLE ?= phi-GRAPE.exe
|
||||||
|
|
||||||
|
default:
|
||||||
|
$(MPICC) $(CPPFLAGS) $(CFLAGS) $(INC) -DDTSCF=$(DTSCF) phi-GRAPE.c -o $(EXECUTABLE) $(LIB)
|
||||||
|
|
||||||
|
yebisu: CPPFLAGS := $(filter-out -DETICS, $(CPPFLAGS))
|
||||||
|
yebisu: default
|
||||||
|
|
||||||
|
clean:
|
||||||
|
rm -f *.o phi-GRAPE.exe
|
||||||
376
act_def_linklist.c
Normal file
376
act_def_linklist.c
Normal file
|
|
@ -0,0 +1,376 @@
|
||||||
|
/**************************************************************
|
||||||
|
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;i<N;i++){
|
||||||
|
|
||||||
|
ttmp = t[i] + dt[i];
|
||||||
|
iii = ind[i];
|
||||||
|
if( m[iii] == 0.0 ) ttmp = t_end + 0.125;
|
||||||
|
|
||||||
|
Tp1 = head;
|
||||||
|
Tp2 = Tp1->NextTNODE;
|
||||||
|
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;i<N;i++)
|
||||||
|
|
||||||
|
CurrT = head->NextTNODE;
|
||||||
|
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;i<n_act;i++){
|
||||||
|
|
||||||
|
iii = ind_act[i];
|
||||||
|
ttmp = t[iii] + dt[iii];
|
||||||
|
|
||||||
|
Tp1 = CurrT;
|
||||||
|
Tp2 = Tp1->NextTNODE;
|
||||||
|
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;i<n_act;i++)
|
||||||
|
|
||||||
|
Tp1 = CurrT;
|
||||||
|
CurrT = CurrT->NextTNODE;
|
||||||
|
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]<tmp) i++;
|
||||||
|
while (tmp<ind_act[j]) j--;
|
||||||
|
|
||||||
|
if (i<=j)
|
||||||
|
{
|
||||||
|
i_swap(&ind_act[i],&ind_act[j]);
|
||||||
|
i++; j--;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
cikl = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
if (l<j) ind_act_sort(l, j);
|
||||||
|
if (i<r) ind_act_sort(i, r);
|
||||||
|
}
|
||||||
|
/**************************************************************/
|
||||||
|
|
||||||
|
/**************************************************************
|
||||||
|
Get active particle list
|
||||||
|
**************************************************************/
|
||||||
|
void get_act_plist()
|
||||||
|
{
|
||||||
|
struct PNODE *Pptr;
|
||||||
|
int i, iii, j,k, itmp, ipt, flag;
|
||||||
|
|
||||||
|
char idcFile[30];
|
||||||
|
FILE *idcomp;
|
||||||
|
|
||||||
|
i=0;
|
||||||
|
Pptr = CurrT->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
|
||||||
|
|
||||||
76
cmd.c
Normal file
76
cmd.c
Normal file
|
|
@ -0,0 +1,76 @@
|
||||||
|
|
||||||
|
int cmd(int kstar, double rstar, double lstar, double Rgal, double *abvmag, double *vmag, double *BV, double *Teff, double *dvmag, double *dBV)
|
||||||
|
{
|
||||||
|
|
||||||
|
double lTeff, BC, kb;
|
||||||
|
double bvc[8], bcc[8];
|
||||||
|
double dbmag;
|
||||||
|
double BCsun, abvmagsun;
|
||||||
|
double rand1, rand2, prand;
|
||||||
|
|
||||||
|
kb = 5.6704E-08*0.5*1.3914E9*0.5*1.3914E9/3.846E26; //Stefan-Boltzmann constant in Lsun Rsun^-2 K^-4
|
||||||
|
|
||||||
|
bvc[0] = -654597.405559323;
|
||||||
|
bvc[1] = 1099118.61158915;
|
||||||
|
bvc[2] = -789665.995692672;
|
||||||
|
bvc[3] = 314714.220932623;
|
||||||
|
bvc[4] = -75148.4728506455;
|
||||||
|
bvc[5] = 10751.803394526;
|
||||||
|
bvc[6] = -853.487897283685;
|
||||||
|
bvc[7] = 28.9988730655392;
|
||||||
|
|
||||||
|
bcc[0] = -4222907.80590972;
|
||||||
|
bcc[1] = 7209333.13326442;
|
||||||
|
bcc[2] = -5267167.04593882;
|
||||||
|
bcc[3] = 2134724.55938336;
|
||||||
|
bcc[4] = -518317.954642773;
|
||||||
|
bcc[5] = 75392.2372207101;
|
||||||
|
bcc[6] = -6082.7301194776;
|
||||||
|
bcc[7] = 209.990478646363;
|
||||||
|
|
||||||
|
BCsun = 0.11; //sun's bolometric correction
|
||||||
|
abvmagsun = 4.83; //sun's absolute V magnitude
|
||||||
|
|
||||||
|
if( rstar && (kstar<14) )
|
||||||
|
{
|
||||||
|
*Teff = pow(lstar/(4.0*Pi*rstar*rstar*kb),0.25);
|
||||||
|
|
||||||
|
if( (*Teff>3000.0) && (*Teff<55000.0) )
|
||||||
|
{
|
||||||
|
lTeff = log10(*Teff);
|
||||||
|
*BV = bvc[0] + bvc[1]*lTeff + bvc[2]*pow(lTeff,2) + bvc[3]*pow(lTeff,3) + bvc[4]*pow(lTeff,4) + bvc[5]*pow(lTeff,5) + bvc[6]*pow(lTeff,6) + bvc[7]*pow(lTeff,7);
|
||||||
|
BC = bcc[0] + bcc[1]*lTeff + bcc[2]*pow(lTeff,2) + bcc[3]*pow(lTeff,3) + bcc[4]*pow(lTeff,4) + bcc[5]*pow(lTeff,5) + bcc[6]*pow(lTeff,6) + bcc[7]*pow(lTeff,7);
|
||||||
|
if(lstar) *abvmag = -2.5*log10(lstar)-BC+BCsun+abvmagsun;
|
||||||
|
*vmag = *abvmag + 5.0*log10(Rgal) - 5.0;
|
||||||
|
|
||||||
|
do{
|
||||||
|
rand1 = 2.0*drand48()-1.0;
|
||||||
|
rand2 = 2.0*drand48()-1.0;
|
||||||
|
} while (rand1*rand1+rand2*rand2 > 1.0);
|
||||||
|
|
||||||
|
prand = sqrt(-2.0*log(rand1*rand1+rand2*rand2)/(rand1*rand1+rand2*rand2));
|
||||||
|
*dvmag = rand1*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2));
|
||||||
|
dbmag = rand2*prand*sqrt(pow(0.02,2) + pow(0.07*pow(10.0, 0.4*(*vmag-25.0)),2));
|
||||||
|
*dBV = *dvmag + dbmag;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
*vmag = 9999.9;
|
||||||
|
*abvmag = 9999.9;
|
||||||
|
*BV = 9999.9;
|
||||||
|
*dvmag = 0.0;
|
||||||
|
*dBV = 0.0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
*Teff = 0.0;
|
||||||
|
*vmag = 9999.9;
|
||||||
|
*abvmag = 9999.9;
|
||||||
|
*BV = 9999.9;
|
||||||
|
*dvmag = 0.0;
|
||||||
|
*dBV = 0.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
return(0);
|
||||||
|
}
|
||||||
143
debug.h
Normal file
143
debug.h
Normal file
|
|
@ -0,0 +1,143 @@
|
||||||
|
/***************************************************************/
|
||||||
|
void d_swap(double *a, double *b)
|
||||||
|
{
|
||||||
|
register double tmp;
|
||||||
|
tmp = *a; *a = *b; *b = tmp;
|
||||||
|
}
|
||||||
|
|
||||||
|
void i_swap(int *a, int *b)
|
||||||
|
{
|
||||||
|
register int tmp;
|
||||||
|
tmp = *a; *a = *b; *b = tmp;
|
||||||
|
}
|
||||||
|
/***************************************************************/
|
||||||
|
|
||||||
|
/***************************************************************/
|
||||||
|
void my_sort(int l, int r, double *arr, int *ind)
|
||||||
|
{
|
||||||
|
|
||||||
|
int i, j, cikl;
|
||||||
|
double tmp;
|
||||||
|
|
||||||
|
i = l; j = r;
|
||||||
|
tmp = arr[(l+r)/2];
|
||||||
|
|
||||||
|
cikl = 1;
|
||||||
|
|
||||||
|
while (cikl)
|
||||||
|
{
|
||||||
|
while (arr[i]<tmp) i++;
|
||||||
|
while (tmp<arr[j]) j--;
|
||||||
|
|
||||||
|
if (i<=j)
|
||||||
|
{
|
||||||
|
d_swap(&arr[i],&arr[j]);
|
||||||
|
i_swap(&ind[i],&ind[j]);
|
||||||
|
i++; j--;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
cikl = 0;
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
if (l<j) my_sort(l, j, arr, ind);
|
||||||
|
if (i<r) my_sort(i, r, arr, ind);
|
||||||
|
}
|
||||||
|
/***************************************************************/
|
||||||
|
|
||||||
|
/***************************************************************/
|
||||||
|
double my_select(int n_1, int n_2, int k, double *arr, int *ind)
|
||||||
|
{
|
||||||
|
|
||||||
|
int i, ir, j, l, mid, a_ind;
|
||||||
|
double a;
|
||||||
|
|
||||||
|
l = n_1;
|
||||||
|
ir = n_2;
|
||||||
|
|
||||||
|
for(;;)
|
||||||
|
{
|
||||||
|
|
||||||
|
if (ir <= l+1)
|
||||||
|
{
|
||||||
|
|
||||||
|
if (ir == l+1 && arr[ir] < arr[l])
|
||||||
|
{
|
||||||
|
d_swap(&arr[l],&arr[ir]);
|
||||||
|
i_swap(&ind[l],&ind[ir]);
|
||||||
|
}
|
||||||
|
|
||||||
|
return(arr[k]);
|
||||||
|
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
|
||||||
|
mid=(l+ir) >> 1;
|
||||||
|
|
||||||
|
d_swap(&arr[mid],&arr[l+1]);
|
||||||
|
i_swap(&ind[mid],&ind[l+1]);
|
||||||
|
|
||||||
|
if (arr[l+1] > arr[ir])
|
||||||
|
{
|
||||||
|
d_swap(&arr[l+1],&arr[ir]);
|
||||||
|
i_swap(&ind[l+1],&ind[ir]);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (arr[l] > arr[ir])
|
||||||
|
{
|
||||||
|
d_swap(&arr[l],&arr[ir]);
|
||||||
|
i_swap(&ind[l],&ind[ir]);
|
||||||
|
}
|
||||||
|
|
||||||
|
if (arr[l+1] > arr[l])
|
||||||
|
{
|
||||||
|
d_swap(&arr[l+1],&arr[l]);
|
||||||
|
i_swap(&ind[l+1],&ind[l]);
|
||||||
|
}
|
||||||
|
|
||||||
|
i = l+1;
|
||||||
|
j = ir;
|
||||||
|
a = arr[l];
|
||||||
|
a_ind = ind[l];
|
||||||
|
|
||||||
|
for (;;)
|
||||||
|
{
|
||||||
|
do i++; while (arr[i] < a);
|
||||||
|
do j--; while (arr[j] > a);
|
||||||
|
if (j < i) break;
|
||||||
|
d_swap(&arr[i],&arr[j]);
|
||||||
|
i_swap(&ind[i],&ind[j]);
|
||||||
|
}
|
||||||
|
|
||||||
|
arr[l] = arr[j];
|
||||||
|
ind[l] = ind[j];
|
||||||
|
arr[j] = a;
|
||||||
|
ind[j] = a_ind;
|
||||||
|
if (j >= k) ir = j-1;
|
||||||
|
if (j <= k) l = i;
|
||||||
|
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
/***************************************************************/
|
||||||
|
|
||||||
|
/***************************************************************/
|
||||||
|
/*
|
||||||
|
double lagrange_radius(double percent)
|
||||||
|
{
|
||||||
|
int Nb;
|
||||||
|
double tmp;
|
||||||
|
|
||||||
|
Nb = (int)(percent * N);
|
||||||
|
|
||||||
|
my_sort(0, N-1, d, ind);
|
||||||
|
|
||||||
|
tmp = my_select(0, N-1, Nb-1, d, ind); d[Nb-1] = tmp;
|
||||||
|
|
||||||
|
return(tmp);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
/***************************************************************/
|
||||||
38
def_DEN.c
Normal file
38
def_DEN.c
Normal file
|
|
@ -0,0 +1,38 @@
|
||||||
|
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
|
||||||
|
void DEF_DEN(double x_gas[][3], double m_gas[], double h_gas[], int sosed_gas[][NB], double DEN_gas[])
|
||||||
|
{
|
||||||
|
|
||||||
|
double Xij, Yij, Zij, Rij, tmp;
|
||||||
|
|
||||||
|
/* SPH style */
|
||||||
|
|
||||||
|
for(i=0;i<N_GAS;i++)
|
||||||
|
{
|
||||||
|
tmp = 0.0;
|
||||||
|
|
||||||
|
for(k=0;k<NB;k++)
|
||||||
|
{
|
||||||
|
j = sosed_gas[i][k];
|
||||||
|
|
||||||
|
Xij = x_gas[i][0]-x_gas[j][0];
|
||||||
|
Yij = x_gas[i][1]-x_gas[j][1];
|
||||||
|
Zij = x_gas[i][2]-x_gas[j][2];
|
||||||
|
|
||||||
|
Rij = sqrt( Xij*Xij + Yij*Yij + Zij*Zij );
|
||||||
|
|
||||||
|
tmp = tmp + 0.5 * m_gas[j] * ( W(Rij,h_gas[i]) + W(Rij,h_gas[j]) );
|
||||||
|
} /* j */
|
||||||
|
|
||||||
|
DEN_gas[i] = tmp;
|
||||||
|
|
||||||
|
} /* i */
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
74
def_DEN_STAR.c
Normal file
74
def_DEN_STAR.c
Normal file
|
|
@ -0,0 +1,74 @@
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
|
||||||
|
void DEF_DEN_STAR(double x_star, double y_star, double z_star, double x_gas[][3], double m_gas[], double h_gas[], double *DEN_star)
|
||||||
|
{
|
||||||
|
|
||||||
|
//int NB_STAR=0;
|
||||||
|
//int sosed_tmp[N_GAS_MAX];
|
||||||
|
|
||||||
|
double Xij, Yij, Zij, Rij, tmp=0.0, tmp2;
|
||||||
|
|
||||||
|
|
||||||
|
/* DEF DEN dlja odnoj STAR chastici */
|
||||||
|
|
||||||
|
//NB_STAR = 0;
|
||||||
|
tmp = 0.0;
|
||||||
|
tmp2 = 0.0;
|
||||||
|
|
||||||
|
for(j=0;j<N_GAS;j++) // GAS
|
||||||
|
{
|
||||||
|
|
||||||
|
tmp2 = SQR(x_star) + SQR(y_star);
|
||||||
|
|
||||||
|
if( (ABS(z_star) < 1e-2) && (tmp2 < 0.2) ) // zamenity na h & (2*R0)^2 = 0.1936
|
||||||
|
{
|
||||||
|
|
||||||
|
Xij = x_star - x_gas[j][0];
|
||||||
|
Yij = y_star - x_gas[j][1];
|
||||||
|
Zij = z_star - x_gas[j][2];
|
||||||
|
|
||||||
|
Rij = sqrt(Xij*Xij + Yij*Yij + Zij*Zij);
|
||||||
|
|
||||||
|
if( Rij < (2.0*h_gas[j]) )
|
||||||
|
{
|
||||||
|
// sosed_tmp[NB_STAR] = j;
|
||||||
|
// NB_STAR++;
|
||||||
|
|
||||||
|
tmp = tmp + m_gas[j] * W(Rij,h_gas[j]);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
// sosed_tmp[j] = j;
|
||||||
|
} /* j */
|
||||||
|
|
||||||
|
//printf("%06d \t %.8E \n",NB_STAR, tmp);
|
||||||
|
|
||||||
|
|
||||||
|
// my_sort(0, N-1, d, sosed_tmp);
|
||||||
|
|
||||||
|
// tmp = my_select(0, N_GAS-1, NB-1, d, sosed_tmp); d[NB-1] = tmp;
|
||||||
|
|
||||||
|
/* SPH style */
|
||||||
|
|
||||||
|
/*
|
||||||
|
tmp = 0.0;
|
||||||
|
|
||||||
|
for(k=0;k<NB_STAR;k++)
|
||||||
|
{
|
||||||
|
j = sosed_tmp[k];
|
||||||
|
|
||||||
|
Rij = sqrt(d[k]);
|
||||||
|
|
||||||
|
tmp = tmp + m[j] * W(Rij,h_gas[j]);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
*DEN_star = tmp;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
85
def_H.c
Normal file
85
def_H.c
Normal file
|
|
@ -0,0 +1,85 @@
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
|
||||||
|
double W(double Rij, double h)
|
||||||
|
{
|
||||||
|
|
||||||
|
double v, v2, v3, tmp=0.0;
|
||||||
|
|
||||||
|
v = Rij/h; v2 = v*v; v3 = v*v*v;
|
||||||
|
|
||||||
|
if( (v >= 0.0) && (v <= 1.0) )
|
||||||
|
{
|
||||||
|
tmp = 1.0 - 1.5*v2 + 0.75*v3;
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
if( (v > 1.0) && (v < 2.0) )
|
||||||
|
{
|
||||||
|
tmp = 0.25*(2.0-v)*(2.0-v)*(2.0-v);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
tmp = 0.0; /* if ( v >= 2.0 ) */
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
tmp = 1.0/(Pi*h*h*h)*tmp;
|
||||||
|
|
||||||
|
return(tmp);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/*************************************************************************/
|
||||||
|
|
||||||
|
void DEF_H(double x_gas[][3], double h_gas[], int sosed_gas[][NB])
|
||||||
|
{
|
||||||
|
|
||||||
|
int sosed_tmp[N_GAS_MAX];
|
||||||
|
double Xij, Yij, Zij;
|
||||||
|
|
||||||
|
for(i=0;i<N_GAS;i++)
|
||||||
|
{
|
||||||
|
|
||||||
|
for(j=0;j<N_GAS;j++)
|
||||||
|
{
|
||||||
|
Xij = x_gas[i][0]-x_gas[j][0];
|
||||||
|
Yij = x_gas[i][1]-x_gas[j][1];
|
||||||
|
Zij = x_gas[i][2]-x_gas[j][2];
|
||||||
|
|
||||||
|
d[j] = Xij*Xij + Yij*Yij + Zij*Zij;
|
||||||
|
|
||||||
|
sosed_tmp[j] = ind_gas[j];
|
||||||
|
} /* j */
|
||||||
|
|
||||||
|
my_sort(0, N_GAS-1, d, sosed_tmp);
|
||||||
|
|
||||||
|
// tmp = my_select(0, N_GAS-1, NB-1, d, sosed_tmp); d[NB-1] = tmp;
|
||||||
|
|
||||||
|
h_gas[i] = sqrt( d[NB-1] )*0.5;
|
||||||
|
|
||||||
|
for(k=0;k<NB;k++) sosed_gas[i][k] = sosed_tmp[k];
|
||||||
|
|
||||||
|
/*
|
||||||
|
h_gas[i] = MAX(h_gas[i], h_min);
|
||||||
|
h_gas[i] = MIN(h_gas[i], h_max);
|
||||||
|
*/
|
||||||
|
|
||||||
|
/*
|
||||||
|
tmp_l = ldiv(i, N/64);
|
||||||
|
|
||||||
|
if(tmp_l.rem == 0)
|
||||||
|
{
|
||||||
|
printf(".");
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
} /* i */
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
|
/*************************************************************************/
|
||||||
185
drag_force.c
Normal file
185
drag_force.c
Normal file
|
|
@ -0,0 +1,185 @@
|
||||||
|
/*****************************************************************************
|
||||||
|
File Name : "Drag Force.c"
|
||||||
|
Contents : Drag force calculus.
|
||||||
|
: Converted & optimized from Chingis Omarov Fortran code.
|
||||||
|
Coded by : Taras Panamarev,
|
||||||
|
: Denis Yurin, Maxim Makukov & Peter Berczik :-)
|
||||||
|
Last redaction : 2012-11-15 15:33
|
||||||
|
*****************************************************************************/
|
||||||
|
|
||||||
|
// main parameters
|
||||||
|
|
||||||
|
const double Md = 0.01;
|
||||||
|
const double s=4, beta_s=0.7, alpha=0.75, h=0.001;
|
||||||
|
const double r_lim2=(0.5*0.5), h_lim2=(5e-3*5e-3);
|
||||||
|
const double gradP_rho = 0.90; // account for -grad(P)/rho in gas disk...
|
||||||
|
|
||||||
|
const double Qtot = 0.01; // 16e3
|
||||||
|
|
||||||
|
/*
|
||||||
|
const double Qtot = 0.00547000; // 8e3
|
||||||
|
const double Qtot = 0.00296975; // 16e3
|
||||||
|
const double Qtot = 0.00193372; // 32e3
|
||||||
|
const double Qtot = 0.00054215; // 128KB
|
||||||
|
*/
|
||||||
|
|
||||||
|
int K;
|
||||||
|
|
||||||
|
double RDISC, //projection of radius vector to the disc plane
|
||||||
|
h_z, // disc profile
|
||||||
|
RDISC2, Z2, RDISC_3over2, R, RR0, RR02, RR0_s, RDISC_5over2, R_zeta, R0_zeta, RRcrit, RRcrit_zeta,Rcrit_zeta, R_zetaR,
|
||||||
|
inner_hole, temp, sqrt_m_bh, pot_drag,
|
||||||
|
sigma0, R02, hR02, Rcrit, R0, zeta,
|
||||||
|
R_DOT, //the first derivative of absolute value of projection
|
||||||
|
COEFDENS,
|
||||||
|
RHO, RHODOT, //disc density and it's firs derivative with respect to time
|
||||||
|
VXDISC, VYDISC,//disc particle velocity components
|
||||||
|
VXRE, VYRE, VZRE, VRE,//the absolute value of relative (star-disc) velocity and it's components
|
||||||
|
coef,
|
||||||
|
FDRAGX, FDRAGY, FDRAGZ,//drag force components
|
||||||
|
DVXRE, DVYRE, DVZRE, DVRE,//first derivative of relative velocity and it's absolute value
|
||||||
|
KFDOT;
|
||||||
|
|
||||||
|
void init_drag_force()
|
||||||
|
{
|
||||||
|
R0 = R_0;
|
||||||
|
Rcrit = R_CRIT;
|
||||||
|
|
||||||
|
sigma0 = (2-alpha)*Md/(TWOPi*sqrt_TWOPi*h*R0);
|
||||||
|
R02 = SQR(R0);
|
||||||
|
hR02 = SQR(h*R0);
|
||||||
|
}
|
||||||
|
|
||||||
|
void calc_drag_force()
|
||||||
|
{
|
||||||
|
|
||||||
|
if (min_t < t_diss_on) return;
|
||||||
|
// if (min_t == t_diss_on) init_drag_force();
|
||||||
|
init_drag_force();
|
||||||
|
|
||||||
|
for (i=0; i<n_act; i++)
|
||||||
|
{
|
||||||
|
|
||||||
|
K = ind_act[i];
|
||||||
|
|
||||||
|
RDISC2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]);
|
||||||
|
Z2 = SQR(x_act_new[i][2]);
|
||||||
|
|
||||||
|
if( (RDISC2 < r_lim2) && (Z2 < h_lim2) )
|
||||||
|
{
|
||||||
|
|
||||||
|
RDISC2 = SQR(x_act_new[i][0]) + SQR(x_act_new[i][1]);
|
||||||
|
R = sqrt(RDISC2 + SQR(x_act_new[i][2]));
|
||||||
|
|
||||||
|
RDISC = sqrt(RDISC2);
|
||||||
|
R_DOT = ( x_act_new[i][0]*v_act_new[i][0] + x_act_new[i][1]*v_act_new[i][1] )/RDISC;
|
||||||
|
|
||||||
|
RR0 = RDISC/R0;
|
||||||
|
RR02 = SQR(RR0);
|
||||||
|
// RR0_s = pow(RR0, s);
|
||||||
|
RR0_s = SQR(RR02);
|
||||||
|
RRcrit = RDISC/Rcrit;
|
||||||
|
|
||||||
|
/*
|
||||||
|
if (RDISC < Rcrit) zeta = 1; else zeta = 0;
|
||||||
|
|
||||||
|
RRcrit_zeta = pow(RRcrit, zeta);
|
||||||
|
R_zeta = pow(RDISC, -2*zeta);
|
||||||
|
R0_zeta = pow(R0, -2*zeta);
|
||||||
|
Rcrit_zeta = pow(Rcrit, -2*zeta);
|
||||||
|
*/
|
||||||
|
|
||||||
|
if (RDISC < Rcrit)
|
||||||
|
{
|
||||||
|
zeta = 1;
|
||||||
|
RRcrit_zeta = RRcrit;
|
||||||
|
R_zeta = 1.0/RDISC2;
|
||||||
|
R0_zeta = 1.0/R02;
|
||||||
|
Rcrit_zeta = 1.0/SQR(Rcrit);
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
zeta = 0;
|
||||||
|
RRcrit_zeta = 1.0;
|
||||||
|
R_zeta = 1.0;
|
||||||
|
R0_zeta = 1.0;
|
||||||
|
Rcrit_zeta = 1.0;
|
||||||
|
}
|
||||||
|
|
||||||
|
R_zetaR = R_zeta/R;
|
||||||
|
|
||||||
|
COEFDENS = pow(RR0, -alpha) * exp(-beta_s*RR0_s)/RRcrit_zeta;
|
||||||
|
RHO = sigma0*COEFDENS*exp(-Z2/(2.0*hR02*SQR(RRcrit_zeta)));
|
||||||
|
|
||||||
|
RHODOT = - R_DOT*(zeta+alpha+beta_s*s*RR0_s )/RDISC
|
||||||
|
- (v_act_new[i][2]*x_act_new[i][2]*R_zeta-zeta*R_zetaR*Z2*R_DOT)/(hR02*Rcrit_zeta);
|
||||||
|
|
||||||
|
RDISC_3over2 = sqrt(RDISC2*RDISC);
|
||||||
|
RDISC_5over2 = RDISC_3over2*RDISC;
|
||||||
|
|
||||||
|
sqrt_m_bh = sqrt(m_bh);
|
||||||
|
|
||||||
|
sqrt_m_bh *= gradP_rho; // account for -grad(P)/rho in gas disk...
|
||||||
|
|
||||||
|
VXDISC = -sqrt_m_bh * x_act_new[i][1]/RDISC_3over2;
|
||||||
|
VYDISC = +sqrt_m_bh * x_act_new[i][0]/RDISC_3over2;
|
||||||
|
|
||||||
|
VXRE = VXDISC - v_act_new[i][0];
|
||||||
|
VYRE = VYDISC - v_act_new[i][1];
|
||||||
|
VZRE = - v_act_new[i][2];
|
||||||
|
|
||||||
|
VRE = sqrt( SQR(VXRE) + SQR(VYRE) + SQR(VZRE) );
|
||||||
|
|
||||||
|
coef = Qtot*RHO*VRE;
|
||||||
|
|
||||||
|
FDRAGX = coef*VXRE;
|
||||||
|
FDRAGY = coef*VYRE;
|
||||||
|
FDRAGZ = coef*VZRE;
|
||||||
|
|
||||||
|
DVXRE = -sqrt_m_bh * v_act_new[i][1]/RDISC_3over2
|
||||||
|
+sqrt_m_bh * 1.5*x_act_new[i][1]*R_DOT/RDISC_5over2
|
||||||
|
-a_act_new[i][0];
|
||||||
|
|
||||||
|
DVYRE = +sqrt_m_bh * v_act_new[i][0]/RDISC_3over2
|
||||||
|
-sqrt_m_bh * 1.5*x_act_new[i][0]*R_DOT/RDISC_5over2
|
||||||
|
-a_act_new[i][1];
|
||||||
|
|
||||||
|
DVZRE = -a_act_new[i][2];
|
||||||
|
|
||||||
|
DVRE = (VXRE*DVXRE + VYRE*DVYRE + VZRE*DVZRE)/VRE;
|
||||||
|
|
||||||
|
KFDOT = Qtot*RHO;
|
||||||
|
|
||||||
|
adot_drag[K][0] = KFDOT*(RHODOT*VXRE*VRE + DVRE*VXRE + DVXRE*VRE);
|
||||||
|
adot_drag[K][1] = KFDOT*(RHODOT*VYRE*VRE + DVRE*VYRE + DVYRE*VRE);
|
||||||
|
adot_drag[K][2] = KFDOT*(RHODOT*VZRE*VRE + DVRE*VZRE + DVZRE*VRE);
|
||||||
|
|
||||||
|
pot_drag = -0.5*( (a_drag[K][0]+FDRAGX)*(x_act_new[i][0]-x_act[i][0]) +
|
||||||
|
(a_drag[K][1]+FDRAGY)*(x_act_new[i][1]-x_act[i][1]) +
|
||||||
|
(a_drag[K][2]+FDRAGZ)*(x_act_new[i][2]-x_act[i][2]) );
|
||||||
|
|
||||||
|
E_sd += m[K]*pot_drag;
|
||||||
|
|
||||||
|
a_drag[K][0] = FDRAGX;
|
||||||
|
a_drag[K][1] = FDRAGY;
|
||||||
|
a_drag[K][2] = FDRAGZ;
|
||||||
|
|
||||||
|
}
|
||||||
|
else
|
||||||
|
{
|
||||||
|
adot_drag[K][0] = 0.0;
|
||||||
|
adot_drag[K][1] = 0.0;
|
||||||
|
adot_drag[K][2] = 0.0;
|
||||||
|
|
||||||
|
a_drag[K][0] = 0.0;
|
||||||
|
a_drag[K][1] = 0.0;
|
||||||
|
a_drag[K][2] = 0.0;
|
||||||
|
} /* if( (RDISC2 < r_lim2) && (Z2 < h_lim2) ) */
|
||||||
|
|
||||||
|
} /* i */
|
||||||
|
|
||||||
|
|
||||||
|
}
|
||||||
|
/***************************************************************/
|
||||||
|
/***************************************************************/
|
||||||
|
/***************************************************************/
|
||||||
76
n_bh.c
Normal file
76
n_bh.c
Normal file
|
|
@ -0,0 +1,76 @@
|
||||||
|
/***************************************************************************/
|
||||||
|
/*
|
||||||
|
Coded by : Peter Berczik
|
||||||
|
Version number : 1.0
|
||||||
|
Last redaction : 2011.II.20. 13:00
|
||||||
|
*/
|
||||||
|
|
||||||
|
int calc_force_n_BH(double m1, double xx1[], double vv1[],
|
||||||
|
double m2, double xx2[], double vv2[],
|
||||||
|
double eps_BH,
|
||||||
|
double *pot_n1, double a_n1[], double adot_n1[],
|
||||||
|
double *pot_n2, double a_n2[], double adot_n2[])
|
||||||
|
{
|
||||||
|
|
||||||
|
/*
|
||||||
|
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
|
||||||
|
|
||||||
|
m2 - mass of the 2 BH
|
||||||
|
xx2[0,1,2] - coordinate of the 2 BH
|
||||||
|
vv2[0,1,2] - velocity of the 2 BH
|
||||||
|
|
||||||
|
eps_BH - force softening, can be even exactly 0.0 !
|
||||||
|
|
||||||
|
OUTPUT
|
||||||
|
|
||||||
|
pot_n1 for the 1 BH
|
||||||
|
a_n1 [0,1,2] for the 1 BH
|
||||||
|
adot_n1 [0,1,2] for the 1 BH
|
||||||
|
|
||||||
|
pot_n2 for the 2 BH
|
||||||
|
a_n2 [0,1,2] for the 2 BH
|
||||||
|
adot_n2 [0,1,2] for the 2 BH
|
||||||
|
*/
|
||||||
|
|
||||||
|
|
||||||
|
int k;
|
||||||
|
|
||||||
|
double r, r2, r3, r4, RP, x[3], v[3];
|
||||||
|
|
||||||
|
|
||||||
|
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]) + SQR(eps_BH);
|
||||||
|
|
||||||
|
r = sqrt(r2);
|
||||||
|
r3 = r2*r;
|
||||||
|
r4 = r2*r2;
|
||||||
|
|
||||||
|
RP = 3.0*(x[0]*v[0]+x[1]*v[1]+x[2]*v[2])/r;
|
||||||
|
|
||||||
|
// Newton pot + acc & jerks
|
||||||
|
|
||||||
|
*pot_n1 = -m2/r;
|
||||||
|
*pot_n2 = -m1/r;
|
||||||
|
|
||||||
|
for(k=0;k<3;k++)
|
||||||
|
{
|
||||||
|
a_n1[k] = -m2*x[k]/r3;
|
||||||
|
a_n2[k] = m1*x[k]/r3;
|
||||||
|
|
||||||
|
adot_n1[k] = -m2*(v[k]/r3 - RP*x[k]/r4);
|
||||||
|
adot_n2[k] = m1*(v[k]/r3 - RP*x[k]/r4);
|
||||||
|
}
|
||||||
|
|
||||||
|
return(0);
|
||||||
|
|
||||||
|
}
|
||||||
|
/***************************************************************************/
|
||||||
7241
phi-GRAPE.c
Normal file
7241
phi-GRAPE.c
Normal file
File diff suppressed because it is too large
Load diff
765
pn_bh.c
Normal file
765
pn_bh.c
Normal file
|
|
@ -0,0 +1,765 @@
|
||||||
|
/***************************************************************************/
|
||||||
|
/*
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
}
|
||||||
|
/***************************************************************************/
|
||||||
765
pn_bh_spin.c
Normal file
765
pn_bh_spin.c
Normal file
|
|
@ -0,0 +1,765 @@
|
||||||
|
/***************************************************************************/
|
||||||
|
/*
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
}
|
||||||
|
/***************************************************************************/
|
||||||
275
star_destr.c
Normal file
275
star_destr.c
Normal file
|
|
@ -0,0 +1,275 @@
|
||||||
|
/*****************************************************************************
|
||||||
|
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<n_end; i++)
|
||||||
|
{
|
||||||
|
|
||||||
|
rsb = sqrt( SQR(x[i][0]-x[i_bh][0]) + SQR(x[i][1]-x[i_bh][1]) + SQR(x[i][2]-x[i_bh][2]) );
|
||||||
|
|
||||||
|
// R_t = (R_s/5.848035) * pow( (2.0*(*m_bh)/m_s), over3 );
|
||||||
|
// if( (r < R_t) && (e_kin < ABS(vir*e_pot)) )
|
||||||
|
|
||||||
|
if( (rsb < R_t) && (m[i] != 0.0) )
|
||||||
|
{
|
||||||
|
|
||||||
|
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
|
||||||
|
out = fopen("accr-data.dat","a");
|
||||||
|
|
||||||
|
fprintf(out,"%04d %04d %04d %04d \t %.6E %.6E %.6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \n",
|
||||||
|
i, ind[i], i_bh, ind[i_bh],
|
||||||
|
time,
|
||||||
|
m[i], *m_bh,
|
||||||
|
x[i][0], x[i][1], x[i][2],
|
||||||
|
x[i_bh][0], x[i_bh][1], x[i_bh][2],
|
||||||
|
v[i][0], v[i][1], v[i][2],
|
||||||
|
v[i_bh][0], v[i_bh][1], v[i_bh][2],
|
||||||
|
pot[i],
|
||||||
|
a[i][0], a[i][1], a[i][2],
|
||||||
|
adot[i][0], adot[i][1], adot[i][2]);
|
||||||
|
|
||||||
|
fclose(out);
|
||||||
|
|
||||||
|
} /* if(myRank == rootRank) */
|
||||||
|
|
||||||
|
|
||||||
|
v_bh[0] = (m[i] * v[i][0] + *m_bh * v[i_bh][0]) / (m[i] + *m_bh);
|
||||||
|
v_bh[1] = (m[i] * v[i][1] + *m_bh * v[i_bh][1]) / (m[i] + *m_bh);
|
||||||
|
v_bh[2] = (m[i] * v[i][2] + *m_bh * v[i_bh][2]) / (m[i] + *m_bh);
|
||||||
|
|
||||||
|
e_kin = m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
|
||||||
|
e_kin += *m_bh * ( SQR(v[i_bh][0]) + SQR(v[i_bh][1]) + SQR(v[i_bh][2]) );
|
||||||
|
e_kin -= (m[i] + *m_bh) * ( SQR(v_bh[0]) + SQR(v_bh[1]) + SQR(v_bh[2]) );
|
||||||
|
e_kin *= 0.5;
|
||||||
|
|
||||||
|
// e_kin = 0.5*m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) - SQR(v[i_bh][0]) - SQR(v[i_bh][1]) - SQR(v[i_bh][2]) );
|
||||||
|
|
||||||
|
e_pot_BH = - *m_bh * m[i] / sqrt( SQR(rsb) + eps_bh2 );
|
||||||
|
|
||||||
|
|
||||||
|
tmp = 0.0;
|
||||||
|
|
||||||
|
for(k=0; k<N_end; k++)
|
||||||
|
{
|
||||||
|
if( k != ind[i] )
|
||||||
|
{
|
||||||
|
|
||||||
|
/* define if k prinadlezhit ACT ili net */
|
||||||
|
|
||||||
|
k_act = 0;
|
||||||
|
for(ii=0; ii<n_end; ii++)
|
||||||
|
{
|
||||||
|
if( k == ind[ii] )
|
||||||
|
{
|
||||||
|
k_act = 1; break;
|
||||||
|
}
|
||||||
|
} /* ii */
|
||||||
|
|
||||||
|
/*
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("0: \t %04d %04d \t %04d %04d \n", i, ind[i], k, k_act);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
if(k_act == 1) /* k prinadlezhit ACT */
|
||||||
|
{
|
||||||
|
|
||||||
|
rkb2 = SQR(x_N[k][0]-x[i_bh][0]) + SQR(x_N[k][1]-x[i_bh][1]) + SQR(x_N[k][2]-x[i_bh][2]);
|
||||||
|
tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
|
||||||
|
|
||||||
|
rks2 = SQR(x_N[k][0]-x[i][0]) + SQR(x_N[k][1]-x[i][1]) + SQR(x_N[k][2]-x[i][2]);
|
||||||
|
tmp += m_N[k]/sqrt( rks2 + eps2 );
|
||||||
|
|
||||||
|
}
|
||||||
|
else /* k ne prinadlezhit ACT */
|
||||||
|
{
|
||||||
|
|
||||||
|
dt_tmp = time - t_N[k];
|
||||||
|
dt2half = over2*dt_tmp*dt_tmp;
|
||||||
|
dt3over6 = over3*dt_tmp*dt2half;
|
||||||
|
|
||||||
|
xp[0] = x_N[k][0] + v_N[k][0]*dt_tmp + a_N[k][0]*dt2half + adot_N[k][0]*dt3over6;
|
||||||
|
xp[1] = x_N[k][1] + v_N[k][1]*dt_tmp + a_N[k][1]*dt2half + adot_N[k][1]*dt3over6;
|
||||||
|
xp[2] = x_N[k][2] + v_N[k][2]*dt_tmp + a_N[k][2]*dt2half + adot_N[k][2]*dt3over6;
|
||||||
|
|
||||||
|
rkb2 = SQR(xp[0]-x[i_bh][0]) + SQR(xp[1]-x[i_bh][1]) + SQR(xp[2]-x[i_bh][2]);
|
||||||
|
tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
|
||||||
|
|
||||||
|
rks2 = SQR(xp[0]-x[i][0]) + SQR(xp[1]-x[i][1]) + SQR(xp[2]-x[i][2]);
|
||||||
|
tmp += m_N[k]/sqrt( rks2 + eps2 );
|
||||||
|
|
||||||
|
} /* if(k_act == 1) */
|
||||||
|
|
||||||
|
/*
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("\t %04d %04d \t %04d \t % .8E \n", k, ind[i], k_act, tmp);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
} /* if( k != ind[i] ) */
|
||||||
|
|
||||||
|
} /* k */
|
||||||
|
|
||||||
|
e_pot = - m[i] * tmp;
|
||||||
|
|
||||||
|
e_corr = e_kin + e_pot_BH + e_pot;
|
||||||
|
|
||||||
|
E_corr += e_corr;
|
||||||
|
|
||||||
|
*m_bh += m[i]; // add the star mass to the BH mass
|
||||||
|
|
||||||
|
*num_bh = *num_bh + 1;
|
||||||
|
|
||||||
|
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("ACCR: %.6E %04d %06d %.4E %.4E %04d %06d %.4E % .4E % .4E % .4E % .4E %.4E %04d \n",
|
||||||
|
time, i_bh, ind[i_bh], rsb, R_t, i, ind[i],
|
||||||
|
e_kin, e_pot_BH, e_pot, e_corr, E_corr, *m_bh, *num_bh);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
v[i_bh][0] = v_bh[0];
|
||||||
|
v[i_bh][1] = v_bh[1];
|
||||||
|
v[i_bh][2] = v_bh[2];
|
||||||
|
|
||||||
|
|
||||||
|
// tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
|
||||||
|
// tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
|
||||||
|
// tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
|
||||||
|
|
||||||
|
|
||||||
|
/* for part "i" */
|
||||||
|
|
||||||
|
x[i][0] = x_max + 1.0*my_rand2();
|
||||||
|
x[i][1] = x_max + 1.0*my_rand2();
|
||||||
|
x[i][2] = x_max + 1.0*my_rand2();
|
||||||
|
|
||||||
|
v[i][0] = v_max*my_rand2();
|
||||||
|
v[i][1] = v_max*my_rand2();
|
||||||
|
v[i][2] = v_max*my_rand2();
|
||||||
|
|
||||||
|
pot[i] = 0.0;
|
||||||
|
|
||||||
|
a[i][0] = a_max*my_rand2();
|
||||||
|
a[i][1] = a_max*my_rand2();
|
||||||
|
a[i][2] = a_max*my_rand2();
|
||||||
|
|
||||||
|
adot[i][0] = adot_max*my_rand2();
|
||||||
|
adot[i][1] = adot_max*my_rand2();
|
||||||
|
adot[i][2] = adot_max*my_rand2();
|
||||||
|
|
||||||
|
m[i] = 0.0;
|
||||||
|
|
||||||
|
t[i] = t[i] + 0.125;
|
||||||
|
|
||||||
|
dt[i] = 0.125;
|
||||||
|
|
||||||
|
} /* if( (r < R_t) ) */
|
||||||
|
|
||||||
|
} /* i */
|
||||||
|
|
||||||
|
|
||||||
|
m[i_bh] = *m_bh;
|
||||||
|
|
||||||
|
/*
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("%06d \t %.6E \n", i_bh, m[i_bh]);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/***************************************************************/
|
||||||
|
/***************************************************************/
|
||||||
|
/***************************************************************/
|
||||||
286
star_destr_ext.c
Normal file
286
star_destr_ext.c
Normal file
|
|
@ -0,0 +1,286 @@
|
||||||
|
/*****************************************************************************
|
||||||
|
File Name : "Star Destr EXT.c"
|
||||||
|
Contents : star "destruction" by tidal field of EXT fixed BH
|
||||||
|
Coded by : Peter Berczik
|
||||||
|
Last redaction : 2010.IX.14 1:43PM
|
||||||
|
*****************************************************************************/
|
||||||
|
|
||||||
|
void star_destr_ext(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 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, eps_bh2, eps2, rsb, rkb2, rks2, xp[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;
|
||||||
|
|
||||||
|
double a_ext_i[3], adot_ext_i[3];
|
||||||
|
|
||||||
|
|
||||||
|
if(time < t_diss_on) return;
|
||||||
|
|
||||||
|
|
||||||
|
eps2 = SQR(eps);
|
||||||
|
|
||||||
|
eps_bh = b_bh;
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
|
||||||
|
n_end = n;
|
||||||
|
N_end = N;
|
||||||
|
|
||||||
|
/*
|
||||||
|
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<n_end; i++)
|
||||||
|
{
|
||||||
|
|
||||||
|
rsb = sqrt( SQR(x[i][0]) + SQR(x[i][1]) + SQR(x[i][2]) );
|
||||||
|
|
||||||
|
// R_t = (R_s/5.848035) * pow( (2.0*(*m_bh)/m_s), over3 );
|
||||||
|
// if( (r < R_t) && (e_kin < ABS(vir*e_pot)) )
|
||||||
|
|
||||||
|
if( (rsb < R_t) )
|
||||||
|
{
|
||||||
|
|
||||||
|
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
|
||||||
|
x_ij = x[i][0];
|
||||||
|
y_ij = x[i][1];
|
||||||
|
z_ij = x[i][2];
|
||||||
|
|
||||||
|
vx_ij = v[i][0];
|
||||||
|
vy_ij = v[i][1];
|
||||||
|
vz_ij = v[i][2];
|
||||||
|
|
||||||
|
r2 = SQR(x_ij) + SQR(y_ij) + SQR(z_ij) + eps_bh2;
|
||||||
|
r = sqrt(r2);
|
||||||
|
|
||||||
|
rv_ij = vx_ij*x_ij + vy_ij*y_ij + vz_ij*z_ij;
|
||||||
|
|
||||||
|
tmp = *m_bh / (r * r2);
|
||||||
|
|
||||||
|
a_ext_i[0] = -tmp * x_ij;
|
||||||
|
a_ext_i[1] = -tmp * y_ij;
|
||||||
|
a_ext_i[2] = -tmp * z_ij;
|
||||||
|
|
||||||
|
adot_ext_i[0] = -tmp * (vx_ij - 3.0*rv_ij * x_ij/r2);
|
||||||
|
adot_ext_i[1] = -tmp * (vy_ij - 3.0*rv_ij * y_ij/r2);
|
||||||
|
adot_ext_i[2] = -tmp * (vz_ij - 3.0*rv_ij * z_ij/r2);
|
||||||
|
|
||||||
|
out = fopen("accr-data.dat","a");
|
||||||
|
|
||||||
|
ii = ind[i];
|
||||||
|
|
||||||
|
#ifndef STARDISK
|
||||||
|
a_drag[ii][0] = 0.0;
|
||||||
|
a_drag[ii][1] = 0.0;
|
||||||
|
a_drag[ii][2] = 0.0;
|
||||||
|
|
||||||
|
adot_drag[ii][0] = 0.0;
|
||||||
|
adot_drag[ii][1] = 0.0;
|
||||||
|
adot_drag[ii][2] = 0.0;
|
||||||
|
#endif // no STARDISK
|
||||||
|
|
||||||
|
fprintf(out, "%04d %04d \t %.6E %.6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \t % .6E % .6E % .6E \n",
|
||||||
|
i, ii,
|
||||||
|
time_cur, *m_bh,
|
||||||
|
x[i][0], x[i][1], x[i][2],
|
||||||
|
v[i][0], v[i][1], v[i][2],
|
||||||
|
a[i][0], a[i][1], a[i][2],
|
||||||
|
adot[i][0], adot[i][1], adot[i][2],
|
||||||
|
a_ext_i[0], a_ext_i[1], a_ext_i[2],
|
||||||
|
adot_ext_i[0], adot_ext_i[1], adot_ext_i[2],
|
||||||
|
a_drag[ii][0], a_drag[ii][1], a_drag[ii][2],
|
||||||
|
adot_drag[ii][0], adot_drag[ii][1], adot_drag[ii][2]);
|
||||||
|
|
||||||
|
fclose(out);
|
||||||
|
|
||||||
|
} /* if(myRank == rootRank) */
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
e_kin = 0.5*m[i] * ( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
|
||||||
|
|
||||||
|
e_pot_BH = - *m_bh * m[i] / sqrt( SQR(rsb) + eps_bh2 );
|
||||||
|
|
||||||
|
|
||||||
|
tmp = 0.0;
|
||||||
|
|
||||||
|
for(k=0; k<N_end; k++)
|
||||||
|
{
|
||||||
|
if( k != ind[i] )
|
||||||
|
{
|
||||||
|
|
||||||
|
/* define if k prinadlezhit ACT ili net */
|
||||||
|
|
||||||
|
k_act = 0;
|
||||||
|
for(ii=0; ii<n_end; ii++)
|
||||||
|
{
|
||||||
|
if(k==ind[ii])
|
||||||
|
{
|
||||||
|
k_act=1; break;
|
||||||
|
}
|
||||||
|
} /* ii */
|
||||||
|
|
||||||
|
/*
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("0: \t %04d %04d \t %04d %04d \n", i, ind[i], k, k_act);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
if(k_act == 1) /* k prinadlezhit ACT */
|
||||||
|
{
|
||||||
|
|
||||||
|
rkb2 = SQR(x_N[k][0]) + SQR(x_N[k][1]) + SQR(x_N[k][2]);
|
||||||
|
tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
|
||||||
|
|
||||||
|
rks2 = SQR(x_N[k][0]-x[i][0]) + SQR(x_N[k][1]-x[i][1]) + SQR(x_N[k][2]-x[i][2]);
|
||||||
|
tmp += m_N[k]/sqrt( rks2 + eps2 );
|
||||||
|
|
||||||
|
}
|
||||||
|
else /* k ne prinadlezhit ACT */
|
||||||
|
{
|
||||||
|
|
||||||
|
dt_tmp = time - t_N[k];
|
||||||
|
dt2half = over2*dt_tmp*dt_tmp;
|
||||||
|
dt3over6 = over3*dt_tmp*dt2half;
|
||||||
|
|
||||||
|
xp[0] = x_N[k][0] + v_N[k][0]*dt_tmp + a_N[k][0]*dt2half + adot_N[k][0]*dt3over6;
|
||||||
|
xp[1] = x_N[k][1] + v_N[k][1]*dt_tmp + a_N[k][1]*dt2half + adot_N[k][1]*dt3over6;
|
||||||
|
xp[2] = x_N[k][2] + v_N[k][2]*dt_tmp + a_N[k][2]*dt2half + adot_N[k][2]*dt3over6;
|
||||||
|
|
||||||
|
rkb2 = SQR(xp[0]) + SQR(xp[1]) + SQR(xp[2]);
|
||||||
|
tmp -= m_N[k]/sqrt( rkb2 + eps_bh2 );
|
||||||
|
|
||||||
|
rks2 = SQR(xp[0]-x[i][0]) + SQR(xp[1]-x[i][1]) + SQR(xp[2]-x[i][2]);
|
||||||
|
tmp += m_N[k]/sqrt( rks2 + eps2 );
|
||||||
|
|
||||||
|
} /* if(k_act == 1) */
|
||||||
|
|
||||||
|
/*
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("\t %04d %04d \t %04d \t % .8E \n", k, ind[i], k_act, tmp);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
*/
|
||||||
|
|
||||||
|
} /* if( k != ind[i] ) */
|
||||||
|
|
||||||
|
} /* k */
|
||||||
|
|
||||||
|
|
||||||
|
e_pot = - m[i] * tmp;
|
||||||
|
|
||||||
|
e_corr = e_kin + e_pot_BH + e_pot;
|
||||||
|
|
||||||
|
E_corr += e_corr;
|
||||||
|
|
||||||
|
*m_bh += m[i]; // add the star mass to the BH mass
|
||||||
|
|
||||||
|
*num_bh = *num_bh + 1;
|
||||||
|
|
||||||
|
|
||||||
|
if(myRank == rootRank)
|
||||||
|
{
|
||||||
|
printf("ACCR: %.4E %.4E %.4E %04d %06d %.4E % .4E % .4E % .4E % .4E %.4E %04d \n",
|
||||||
|
time, rsb, R_t, i, ind[i], e_kin, e_pot_BH, e_pot, e_corr, E_corr, *m_bh, *num_bh);
|
||||||
|
fflush(stdout);
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
// tmp_v = sqrt( SQR(v[i][0]) + SQR(v[i][1]) + SQR(v[i][2]) );
|
||||||
|
// tmp_a = sqrt( SQR(a[i][0]) + SQR(a[i][1]) + SQR(a[i][2]) );
|
||||||
|
// tmp_adot = sqrt( SQR(adot[i][0]) + SQR(adot[i][1]) + SQR(adot[i][2]) );
|
||||||
|
|
||||||
|
|
||||||
|
/* for part "i" */
|
||||||
|
|
||||||
|
x[i][0] = x_max + 1.0*my_rand2();
|
||||||
|
x[i][1] = x_max + 1.0*my_rand2();
|
||||||
|
x[i][2] = x_max + 1.0*my_rand2();
|
||||||
|
|
||||||
|
v[i][0] = v_max*my_rand2();
|
||||||
|
v[i][1] = v_max*my_rand2();
|
||||||
|
v[i][2] = v_max*my_rand2();
|
||||||
|
|
||||||
|
pot[i] = 0.0;
|
||||||
|
|
||||||
|
a[i][0] = a_max*my_rand2();
|
||||||
|
a[i][1] = a_max*my_rand2();
|
||||||
|
a[i][2] = a_max*my_rand2();
|
||||||
|
|
||||||
|
adot[i][0] = adot_max*my_rand2();
|
||||||
|
adot[i][1] = adot_max*my_rand2();
|
||||||
|
adot[i][2] = adot_max*my_rand2();
|
||||||
|
|
||||||
|
m[i] = 0.0;
|
||||||
|
|
||||||
|
t[i] = t[i] + 0.125;
|
||||||
|
|
||||||
|
dt[i] = 0.125;
|
||||||
|
|
||||||
|
} /* if( (r < R_t) ) */
|
||||||
|
|
||||||
|
} /* i */
|
||||||
|
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
/***************************************************************/
|
||||||
|
/***************************************************************/
|
||||||
|
/***************************************************************/
|
||||||
Loading…
Add table
Add a link
Reference in a new issue