Created
September 4, 2014 11:39
-
-
Save cmaureir/64d8e3505ef6dccaee11 to your computer and use it in GitHub Desktop.
Accretion
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| // ... | |
| if(P[j].Mass > 0) | |
| { | |
| if(mass > 0) | |
| { | |
| dx = pos[0] - P[j].Pos[0]; | |
| dy = pos[1] - P[j].Pos[1]; | |
| dz = pos[2] - P[j].Pos[2]; | |
| r2 = dx * dx + dy * dy + dz * dz; | |
| //if(r2 < h_i2) | |
| // Inside the softening | |
| if(r2 < 100 * All.SofteningTable[5] * All.SofteningTable[5]) | |
| { | |
| if(P[j].Type == 0) | |
| { | |
| /* here we have a gas particle */ | |
| r = sqrt(r2); | |
| hinv = 1 / h_i; | |
| #ifndef TWODIMS | |
| hinv3 = hinv * hinv * hinv; | |
| #else | |
| hinv3 = hinv * hinv / boxSize_Z; | |
| #endif | |
| u = r * hinv; | |
| if(u < 0.5) | |
| wk = hinv3 * (KERNEL_COEFF_1 + KERNEL_COEFF_2 * (u - 1) * u * u); | |
| else | |
| wk = hinv3 * KERNEL_COEFF_5 * (1.0 - u) * (1.0 - u) * (1.0 - u); | |
| #ifdef SWALLOWGAS | |
| #ifdef SWALLOW_BH_BINARY | |
| // Omitting primary | |
| if (id != 200001) | |
| { | |
| // BHvel - GASvel | |
| double vdx = velocity[0] - P[j].Vel[0]; | |
| double vdy = velocity[1] - P[j].Vel[1]; | |
| double vdz = velocity[2] - P[j].Vel[2]; | |
| double v2 = vdx*vdx + vdy*vdy + vdz*vdz; | |
| double normrel = v2 + csnd*csnd; | |
| double r_bound = (2 * All.G * mass) / normrel; | |
| double r_fix = 0.05; | |
| printf("====> (%d, %d) |ID:%d (%d)| r=%.5e r_fix=%.5e r_bound=%.5e mass=%.5e norm=%.5e\n", | |
| target, mode, | |
| P[target].ID, P[target].Type, | |
| r, r_fix, r_bound, | |
| mass, normrel); | |
| // Recipe from Amaro-Seoane/Dotti/Colpi paper. | |
| #ifdef RFIX | |
| if (r <= 0.5 * r_fix) | |
| { | |
| P[j].SwallowID = id; | |
| printf("[0] Particle %d will be swallow by %d\n", P[j].ID, P[j].SwallowID); | |
| // Redefining the epsilon to have 2/3 of the material | |
| // accreted. | |
| epsilon = 1.0/3.0; | |
| } | |
| else if ( 0.5 * r_fix < r && r <= r_fix) | |
| { | |
| P[j].SwallowID = id; | |
| printf("[1] Particle %d will be swallow by %d\n", P[j].ID, P[j].SwallowID); | |
| // Redefining the epsilon to have 1/3 of the material | |
| // accreted. | |
| epsilon = 2.0/3.0; | |
| } | |
| #else | |
| if (r < alpha * r_bound) | |
| { | |
| P[j].SwallowID = id; | |
| printf("[3] Particle %d will be swallow by %d\n", P[j].ID, P[j].SwallowID); | |
| } | |
| #endif | |
| } | |
| #endif // SWALLOW_BH_BINARY | |
| // Old stochastic accretion | |
| // double p = 0.0, w = 0.0; | |
| //#ifndef LT_BH_ACCRETE_SLICES | |
| // if((bh_mass - mass) > 0 && rho > 0) | |
| // p = (bh_mass - mass) * wk / rho; | |
| //#else | |
| // if((bh_mass - mass) > 0 && rho > 0) | |
| // p = (bh_mass - mass) * wk / rho * (All.NBHslices - SphP[j].NSlicesSwallowed); | |
| //#endif | |
| // else | |
| // p = 0; | |
| // | |
| // w = get_random_number(P[j].ID); | |
| // if(w < p) | |
| // { | |
| // if(P[j].SwallowID < id) | |
| // { | |
| // P[j].SwallowID = id; | |
| // } | |
| // } | |
| // | |
| #endif |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment