Nbody6
dgcore.f
Go to the documentation of this file.
1 ***
2  SUBROUTINE dgcore(kw1,kw2,kw3,m1,m2,m3,ebinde)
3 *
4 * A routine to determine the outcome of a collision or coalescence
5 * of two degenerate cores.
6 * Entered with kw1,kw2 = 2 or 3 with M <= Mflash, 6, 10, 11 or 12
7 *
8  implicit none
9 *
10  integer kw1,kw2,kw3
11 *
12  real*8 m1,m2,m3,ebinde
13  real*8 r1,r2,r3,mhe,mc,mne,ebindi,ebindf,deleb,de,enuc
14  real*8 temp,x,y,m0,mflash
15  real*8 cvhe,cvc,cvne
16  parameter(cvhe=3.1d+07,cvc=8.27d+06,cvne=7.44d+06)
17  real*8 ehe,ec,ene
18  parameter(ehe=5.812d+17,ec=2.21d+17,ene=2.06d+17)
19  real*8 the,tc,gmr,mch
20  parameter(the=1.0d+08,tc=1.0d+09,gmr=1.906d+15,mch=1.44d0)
21 *
22  real*8 corerd
23  external corerd
24 *
25 * Calculate the core radii setting m0 < mflash using dummy values as we
26 * know it to be true if kw = 2 or 3.
27  m0 = 1.d0
28  mflash = 2.d0
29  r1 = corerd(kw1,m1,m0,mflash)
30  r2 = corerd(kw2,m2,m0,mflash)
31  r3 = corerd(kw3,m3,m0,mflash)
32 * Calculate the initial binding energy of the two seperate cores and the
33 * difference between this and the final binding energy of the new core.
34  ebindi = m1*m1/r1 + m2*m2/r2
35  ebindf = m3*m3/r3
36  deleb = abs(ebindi - ebindf)
37 * If an envelope is present reduce its binding energy by the amount
38 * of energy liberated by the coalescence.
39  ebinde = max(0.d0,ebinde - deleb)
40  if(kw1.gt.3) goto 90
41 * Distribute the mass into core mass groups where mhe represents Helium
42 * core mass, mc represents Carbon core mass and mne represents a Neon
43 * core mass or any mass that is all converted Carbon.
44  mhe = 0.d0
45  mc = 0.d0
46  mne = 0.d0
47  if(kw1.le.3.or.kw1.eq.10)then
48  mhe = mhe + m1
49  elseif(kw1.eq.12)then
50  mne = mne + m1
51  else
52  mc = mc + m1
53  endif
54  if(kw2.le.3.or.kw2.eq.10)then
55  mhe = mhe + m2
56  elseif(kw2.eq.12)then
57  mne = mne + m2
58  else
59  mc = mc + m2
60  endif
61 * Calculate the temperature generated by the merging of the cores.
62  temp = (deleb/(cvhe*mhe+cvc*mc+cvne*mne))*gmr
63 *
64 * To decide if He is converted to C we use:
65 * 3He4 -> C12 , T > 10^8 K , 7.274 Mev released,
66 * to decide if C is converted further we use:
67 * 2C12 -> Ne20 + alpha , T > 10^9 K , 4.616 Mev released.
68 * and to decide if O is converted further we use:
69 * 2O16 -> P31 + p , T > 10^9 K , 7.677 Mev released.
70 * To obtain the heat capacity of an O/Ne WD and to gain an idea of the
71 * energy released from the further processing of an O/Ne WD we use:
72 * 2Ne20 + gamma -> O16 + Mg24 +gamma , T > 10^9 K , 4.583 Mev released.
73 * For a CO core the composition is assumed to be 20% C, 80% O and for
74 * an ONe core 80% O, 20% Ne.
75 *
76 * Decide if conversion of elements can take place.
77 * if(temp.gt.the)then
78  x = 1.d0
79 * else
80 * x = 0.d0
81 * endif
82 * if(temp.gt.tc)then
83 * y = 1.d0
84 * else
85  y = 0.d0
86 * endif
87 * Calculate the nuclear energy generated from any element conversion.
88  enuc = (x*ehe*mhe + y*(ec*mc + ene*mne))/gmr
89 * Calculate the difference between the binding energy of the star
90 * (core + envelope) and the nuclear energy. The new star will be
91 * destroyed in a SN if dE < 0.
92  de = (ebindf + ebinde) - enuc
93 * If the star survives and an envelope is present then reduce the
94 * envelope binding energy by the amount of liberated nuclear energy.
95 * The envelope will not survive if its binding energy is reduced to <= 0.
96  if(de.ge.0.d0) ebinde = max(0.d0,ebinde - enuc)
97 * Now determine the final evolution state of the merged core.
98  if(de.lt.0.d0) kw3 = 15
99  if(kw3.eq.3)then
100  if(x.gt.0.5d0)then
101  kw3 = 6
102  elseif(ebinde.le.0.d0)then
103  kw3 = 10
104  endif
105  elseif(kw3.eq.4)then
106  if(x.gt.0.5d0)then
107  kw3 = 6
108  elseif(ebinde.le.0.d0)then
109  kw3 = 7
110  endif
111  endif
112  if(kw3.eq.6.and.y.lt.0.5d0)then
113  if(ebinde.le.0.d0) kw3 = 11
114  elseif(kw3.eq.6.and.y.gt.0.5d0)then
115  if(ebinde.le.0.d0) kw3 = 12
116  endif
117  if(kw3.eq.10.and.x.gt.0.5d0) kw3 = 11
118  if(kw3.eq.11.and.y.gt.0.5d0) kw3 = 12
119  if(kw3.ge.10.and.kw3.le.12.and.m3.ge.mch) kw3 = 15
120 *
121  if(kw3.eq.15) m3 = 0.d0
122  90 continue
123 *
124  return
125  end
126 ***