c************************************************************************** SUBROUTINE ppcno3an(t,ro,comp,dcomp,jac,deriv,fait, 1 epsilon,et,ero,ex,hhe,be7e,b8e,n13e,o15e,f17e) c routine PRIVATE du module mod_nuc c He burning: cycles PP, CNO et 3 alpha, avec neutrons, extension de ppcno3a9 c cf. Clayton p. 380, 392 et 430 + cours B.Pichon c éléments pris en compte: c H1, He3, He4, C12, C13, N14, N15, O16, O17, O18, Ne20, Ne22, Ne23, Na23, Mg24, c Mg25, Mg26, Al27, Si28, Si29, Ex c Ex est l'élément fictif complément, il n'intéresse que la diffusion c H2, Li7, Be7 à l'équilibre c Nombre d'espèces 12, nombre d'isotopes 22 c un premier appel a rq_reac initialise et definit le nb. c d'éléments chimiques pour lesquels les reac. nuc. sont tabulees c dans ppcno93a9 on ajoute Ex, soit nchim+1 c on utilise une table ppcno+3a à 9 éléments créée par tab_reac c dont on doit indiquer le nom dans le fichier de données *.don c Auteur: P. Morel, Département J.D. Cassini, O.C.A., CESAM2k c entrées : c t : température cgs c ro : densité cgs c comp : abondances par mole c deriv=.true. : on calcule le jacobien c fait=1 : initialisation de la composition chimique c =2 : calcul de dcomp et jacobien si deriv c =3 : énergie nucléaire et dérivées / t et ro c =4 : production de neutrinos c sorties c dcomp : dérivée temporelle (unité de temps : 10**6 ans) c jac : jacobien (unité de temps : 10**6 ans) c epsilon, et, ero, ex : énergie thermonucléaire (unité de temps : s) c : et dérivées /t, ro ,X c hhe, be7e, b8e, n13e, o15e, f17e : nombre de neutrinos g/s c hhe réaction : H1(p,e+ nu)H2 c be7e réaction : Be7(e-,nu g)Li7 c b8e réaction : B8(,e+ nu)Be8 c n13e réaction : N13(,e+ nu)C13 c o15e réaction : O15(e+,nu)N15 c f17e réaction : F17(,e+ nu)O17 c réactions thermonucléaires de la combustion de l'hélium, cours de B.Pichon p.323 c r(1) : réaction H1(p,e+ nu)H2 (1) PP c r(2) : réaction H2(p,g)He3 (2) c r(3) : réaction He3(He3,2p)He4 (3) c r(4) : réaction He3(a,g)Be7 (4) c r(5) : réaction Li7(p,a)He4 (5) c r(6) : réaction Be7(e-,nu g)Li7 (6) c r(7) : réaction Be7(p,g)B8(e+ nu)Be8(a)He4 (7) c r(8) : réaction C12(p,g)N13(,e+ nu)C13 (8) CNO c r(9) : réaction C13(p,g)N14 (9) c r(10) : réaction N14(p,g)O15(e+,nu)N15 (10) c r(11) : réaction N15(p,g)O16 (11) c r(12) : réaction N15(p,a)C12 (12) c r(13) : réaction O16(p,g)F17(,e+ nu)O17 (13) c r(14) : réaction O17(p,a)N14 (14) c r(15) : réaction He4(2a,g)C12 (15) 3 alpha c r(16) : réaction C12(a,g)O16 (16) c r(17) : réaction O16(a,g)Ne20 (17) c r(18) : N14(a,g)F18(e+ nu)O18 (32) c r(19) : O18(a,g)Ne22 (33) c r(20) : Ne22(a,g)Mg26 (35) c r(21) : Ne22(a,n)Mg25 (34) c r(22) : Ne22(n,g)Ne23(e- nu)Na23 (100) c r(23) : Na23(n,g)Na24(e- nu)Mg24 neutron poison (75) c r(24) : Ne20(n,g)Ne21 neutron poison (72) c r(25) : Ne21(a,n)Mg24 (37) c r(26) : Mg24(n,g)Mg25 neutron poison (76) c r(27) : Mg25(n,g)Mg26 neutron poison (77) c r(28) : Mg26(n,g)Mg27(e- nu)Al27 neutron poison (101) c r(29) : Al27(n,g)Al28(e- nu)Si28 neutron poison (81) c r(30) : Si28(n,g)Si29 neutron poison (84) c r(31) : C13(a,n)O16 (30) c r(32) : Mg25(a,n)Si28 (78) c r(33) : Mg26(a,n)Si29 (74) c r(34) : O18(p,a)N15 (28) c r(35) : réaction O17(p,g)F18(e+ nu)O18 (24) c r(36) : réaction Al27(p,a)Mg24 (58) c r(37) : réaction Al27(p,g)Si28 (57) c r(38) : réaction Na23(p,g)Mg24 (22) c r(39) : réaction Ne22(p,g)Na23 (63) c r(40) : réaction Ne20(g,a)O16 (47) c r(41) : réaction Mg24(g,a)Ne20 (61) c r(42) : réaction Ne20(a,g)Mg24 (29) c r(43) : réaction Mg24(a,g)Si28 (56) c r(44) : réaction C12(C12,g)Mg24 (18) c r(45) : réaction C12(C12,p)Na23 (20) c r(46) : réaction C12(C12,a)Ne20 (21) c r(47) : réaction O16(C12,g)Si28 (48) c r(48) : réaction O16(C12,p)Al27 (50) c r(49) : réaction O16(C12,a)Mg24 (51) c r(50) : réaction O16(O16,a)Si28 (55) c r(51) : réaction C12(C12,n)Mg23(e+ nu)Na23 (19) c réaction O16(O16,g)S32 (52) c réaction O16(O16,n)S31 (53) c réaction O16(O16,p)P31 (54) c indices des éléments c H1 : 1 c He3 : 2 c He4 : 3 i_he4 c C12 : 4 c C13 : 5 c N14 : 6 c N15 : 7 c O16 : 8 c O17 : 9 c O18 : 10 c Ne20 : 11 c Ne21 : 12 c Ne22 : 13 c Na23 : 14 c Mg24 : 15 c Mg25 : 16 c Mg26 : 17 c Al27 : 18 c Si28 : 19 c Si29 : 20 c n : 21 c Ex : 22 i_ex c---------------------------------------------------------------------- USE mod_donnees, ONLY : ab_ini, ah, amu, ihe4, i_ex, langue, nchim, 1 nom_elem, nom_xheavy, nucleo, secon6, t_inf, x0, y0, zi, z0 USE mod_kind USE mod_numerique, ONLY : gauss_band IMPLICIT NONE INTEGER, INTENT(in) :: fait LOGICAL, INTENT(in) :: deriv REAL (kind=dp), INTENT(in):: t, ro REAL (kind=dp), INTENT(inout), DIMENSION(:) :: comp REAL (kind=dp), INTENT(out), DIMENSION(:,:) :: jac REAL (kind=dp), INTENT(out), DIMENSION(:) :: dcomp, ex, epsilon REAL (kind=dp), INTENT(out) :: et, ero, hhe, be7e, b8e, n13e, 1 o15e, f17e REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: drx, dqx REAL (kind=dp), ALLOCATABLE, DIMENSION(:,:) :: a, b REAL (kind=dp), ALLOCATABLE, SAVE, DIMENSION(:) :: anuc, comp_dex, 1 dmuex, dh2x, denx, dbe7x, dli7x, drt, dro, r, q, dqt, dqo REAL (kind=dp) :: mue, nbz, den, be7, h2, li7, dh2t, dh2ro, 1 dent, denro, dbe7t, dbe7ro, dli7t, dli7ro, 2 mass_ex, charge_ex, sum_a INTEGER, ALLOCATABLE, DIMENSION(:) :: indpc INTEGER :: i, j LOGICAL :: inversible CHARACTER (len=2) :: text c-------------------------------------------------------------------------- 2000 FORMAT(8es10.3) 2001 FORMAT(5es15.8) 2002 FORMAT(11es8.1) c initialisations SELECT CASE(fait) CASE(0) c définition de nchim: nombre d'éléments chimiques dont on calcule l'abondance nchim=21+1 c appel d'initialisation pour tabulation des réactions nucléaires c allocations fictives ALLOCATE(drx(1,1),dqx(1,1),r(1),drt(1),dro(1),q(1), 1 dqt(1),dqo(1),dmuex(1)) CALL rq_reac(comp,1.d7,1.d0,r,drt,dro,drx,q,dqt,dqo,dqx,mue,dmuex) DEALLOCATE(dqx,drx) ; ALLOCATE(dqx(nreac,nchim),drx(nreac,nchim)) c--------------------------------------------------------------------- CASE(1) c détermination des abondances initiales selon les données de ab_ini c He3+He4=Y0 c Z0 = somme des éléments plus lourds que hélium c dans Z rapports en nombre CALL abon_ini c Ex : élément fictif moyenne des éléments # CNO charge_ex=0.d0 ; mass_ex=0.d0 ; sum_a=0.d0 DO i=5,nelem_ini !à partir de B=5 sauf les espèces explicitement utilisés IF(elem(i) == ' C')CYCLE IF(elem(i) == ' N')CYCLE IF(elem(i) == ' O')CYCLE IF(elem(i) == 'Ne')CYCLE IF(elem(i) == 'Na')CYCLE IF(elem(i) == 'Mg')CYCLE IF(elem(i) == 'Al')CYCLE IF(elem(i) == 'Si')CYCLE IF(elem(i) == 'n')CYCLE charge_ex=charge_ex+c(i)*ab(i); mass_ex=mass_ex+m(i)*ab(i) sum_a=sum_a+ab(i) ENDDO charge_ex=NINT(charge_ex/sum_a) ; mass_ex=NINT(mass_ex/sum_a) WRITE(text,10)NINT(mass_ex) 10 FORMAT(i2) c élément fictif Ex, indice 21 i_ex=22 !indice de l'élément chimique reliquat nucleo(i_ex)=mass_ex !nucleo de l'élément chimique reliquat zi(i_ex)=charge_ex !charge de l'élément chimique reliquat i=NINT(charge_ex) nom_elem(i_ex)=elem(i)//text !nom elem. chim. rel. nom_xheavy=nom_elem(i_ex) SELECT CASE(langue) CASE('english') WRITE(*,1023)NINT(mass_ex),NINT(charge_ex),TRIM(nom_elem(i_ex)) WRITE(2,1023)NINT(mass_ex),NINT(charge_ex),TRIM(nom_elem(i_ex)) 1023 FORMAT('Fictitious species of mass : ',i3,'and charge :',i3,' : ',a) CASE DEFAULT WRITE(*,23)NINT(mass_ex),NINT(charge_ex),TRIM(nom_elem(i_ex)) WRITE(2,23)NINT(mass_ex),NINT(charge_ex),TRIM(nom_elem(i_ex)) 23 FORMAT('Elément fictif de masse : ',i3,' et de charge :',i3,' : ',a) END SELECT c PRINT*,nchim ; WRITE(*,2000)nucleo(1:nchim) c PRINT*,'abon_rela' ; WRITE(*,2000)abon_rela c PRINT*,nom_elem ; PAUSE'abon_rela' c détermination des abondances initiales, a(équation,isotope) c cf. Eq. (7,97) p. 135 de la NOTICE ALLOCATE(a(nchim,nchim),indpc(nchim),b(1,nchim)) a=0.d0 ; b=0.d0 ; indpc=1 c x1*nu1=X=H1 a(1,1)=nucleo(1) !H1 b(1,1)=x0 c x2*nu2+x3*nu3=Y=He4 a(2,2)=nucleo(2) !He3 a(2,3)=nucleo(3) !He4 b(1,2)=y0 c abon_rela : abondance relative des métaux dans Z, défini dans abon_ini DO j=4,nchim c SUM(j > He4)xi*nui=Z a(3,j)=nucleo(j) c coeff. de xj dans -Xi/Z SUM(j > He4)xj a(4,j)=-abon_rela(6) !somme comp(i) C, C/Z a(5,j)=-abon_rela(7) !somme comp(i) N, N/Z a(6,j)=-abon_rela(8) !somme comp(i) O, O/Z a(7,j)=-abon_rela(10) !somme comp(i) Ne, Ne/Z a(8,j)=-abon_rela(11) !somme comp(i) Na, Na/Z a(9,j)=-abon_rela(12) !somme comp(i) Mg, Mg/Z a(10,j)=-abon_rela(13)!somme comp(i) Al, Al/Z a(11,j)=-abon_rela(14)!somme comp(i) Si, Si/Z a(12,j)=-abon_rela(0) !somme comp(i) n, n/Z ENDDO b(1,3)=z0 !Z c coeff. de x_k.. dans [x_k+x_k+1..] -Xi/Z SUM(j > He4)xj = 0 a(4,4)=a(4,4)+1.d0 !C12 a(4,5)=a(4,5)+1.d0 !C13 a(5,6)=a(5,6)+1.d0 !N14 a(5,7)=a(5,7)+1.d0 !N15 a(6,8)=a(6,8)+1.d0 !O16 a(6,9)=a(6,9)+1.d0 !O17 a(6,10)=a(6,10)+1.d0 !O18 a(7,11)=a(7,11)+1.d0 !Ne20 a(7,12)=a(7,12)+1.d0 !Ne21 a(7,13)=a(7,13)+1.d0 !Ne22 a(8,14)=a(8,14)+1.d0 !Na23 a(9,15)=a(9,15)+1.d0 !Mg24 a(9,16)=a(9,16)+1.d0 !Mg25 a(9,17)=a(9,17)+1.d0 !Mg26 a(10,18)=a(10,18)+1.d0!Al27 a(11,19)=a(11,19)+1.d0!Si28 a(11,20)=a(11,20)+1.d0!Si29 a(12,21)=a(12,21)+1.d0!n c équations des rapports isotopiques a(13,2)=1.d0 !He3 a(13,3)=-he3she4z !He3/He4 avec H2 dans He3 a(14,5)=1.d0 !C13 a(14,4)=-c13sc12 !C13/C12 a(15,7)=1.d0 !N15 a(15,6)=-n15sn14 !N15/N14 a(16,9)=1.d0 !O17 a(16,8)=-o17so16 !O17/O16 a(17,10)=1.d0 !O18 a(17,8)=-o18so16 !O18/O16 a(18,12)=1.d0 !Ne21 a(18,11)=-ne21sne20 !Ne21/Ne20 a(19,13)=1.d0 !Ne22 a(19,11)=-ne22sne20 !Ne22/Ne20 a(20,16)=1.d0 !Mg25 a(20,15)=-mg25smg24 !Mg25/Mg24 a(21,17)=1.d0 !Mg26 a(21,16)=-mg26smg25 !Mg26/Mg25 a(22,20)=1.d0 !Si29 a(22,19)=-si29ssi28 !Si29/Si28 c PRINT*,nchim c DO i=1,nchim c WRITE(*,2002)a(i,1:nchim),b(1,i) c ENDDO CALL gauss_band(a,b,indpc,nchim,nchim,nchim,1,inversible) IF(.NOT.inversible)THEN PRINT*,'ppcno3an, matrice du calcul des abon. non inversible' PRINT*,'ARRET' ; STOP ENDIF c allocations diverses DEALLOCATE(drt,dro,r,q,dqt,dqo,dmuex) ALLOCATE(ab_ini(nchim),anuc(nchim), 1 comp_dex(nchim),dbe7x(nchim),denx(nchim),dh2x(nchim), 2 dli7x(nchim),dmuex(nchim),dqo(nreac),dqt(nreac),dro(nreac), 3 drt(nreac),q(nreac),r(nreac)) c abondances initiales et abondances négligeables comp(:)=b(1,:) ab_ini=comp*nucleo c WRITE(*,2000)comp ; PAUSE'comp ppcno3an' c nombre/volume des métaux dans Z nbz=SUM(comp(ihe4+1:nchim)) c abondances en DeX, H=12 comp_dex=12.d0+LOG10(comp/comp(1)) c écritures WRITE(2,2) ; WRITE(*,2) 2 FORMAT(/,'Réac. thermo. de la combustion de He, cours B.Pichon p.323',/) WRITE(2,3)nreac ; WRITE(*,3)nreac 3 FORMAT('nombre de réactions : ',i3) WRITE(2,4)nreac ; WRITE(*,4)nchim 4 FORMAT('nombre d''éléments chimiques : ',i3) WRITE(2,20)x0,y0,z0,z0/x0 ; WRITE(*,20)x0,y0,z0,z0/x0 20 FORMAT(/,'abondances initiales déduites de:',/, 1 'X0=',es10.3,', Y0=',es10.3,', Z0=',es10.3,/,'Z0/X0=',es10.3, 2 ', H1=X0, He3+He4=Y0',/, 3'Z0=C12+C13+N14+N15+O16+O17+O18+Ne20+Ne21Ne22+Na23+Mg24+Mg25+Mg26+Al27+Si28+Si29+n+Ex') WRITE(2,1)ab_ini ; WRITE(*,1)ab_ini 1 FORMAT(/,'abondances initiales/gramme:',/, 1 'H1:',es10.3,', He3:',es10.3,', He4:',es10.3,', C12:',es10.3, 2 ', C13:',es10.3,/'N14:',es10.3,', N15:',es10.3,', O16:',es10.3, 3 ', O17:',es10.3,', O18:',es10.3,/,'Ne20:',es10.3,', Ne21:',es10.3, 4 ', Ne22:',es10.3,', Na23:',es10.3,', Mg24:',es10.3,/,'Mg25:',es10.3, 5 ', Mg26:',es10.3,', Al27:',es10.3,', Si28:',es10.3,', Si29:',es10.3,/, 6 'n:',es10.3,', Ex:',es10.3) WRITE(2,9)comp_dex ; WRITE(*,9)comp_dex 9 FORMAT(/,'Abondances initiales en nombre: 12+Log10(Ni/Nh)',/, 1 'H1:',es10.3,', He3:',es10.3,', He4:',es10.3,', C12:',es10.3, 2 ', C13:',es10.3,/'N14:',es10.3,', N15:',es10.3,', O16:',es10.3, 3 ', O17:',es10.3,', O18:',es10.3,/,'Ne20:',es10.3,', Ne21:',es10.3, 4 ', Ne22:',es10.3,', Na23:',es10.3,', Mg24:',es10.3,/,'Mg25:',es10.3, 5 ', Mg26:',es10.3,', Al27:',es10.3,', Si28:',es10.3,', Si29:',es10.3,/, 6 'n:',es10.3,', Ex:',es10.3) WRITE(2,21)(comp(4)+comp(5))/nbz, !C/Z 1 (comp(6)+comp(7))/nbz, !N/Z 2 (comp(8)+comp(9)+comp(10))/nbz, !O/Z 3 (comp(11)+comp(12))/nbz, !Ne/Z 4 comp(13)/nbz, !Na/Z 5 (comp(14)+comp(15)+comp(16))/nbz, !Mg/Z 6 comp(17)/nbz, !Al/Z 7 (comp(18)+comp(19))/nbz, !Si/Z 8 comp(20)/nbz, !n/Z 9 comp(21)/nbz !Ex/Z WRITE(*,21)(comp(4)+comp(5))/nbz, !C/Z 1 (comp(6)+comp(7))/nbz, !N/Z 2 (comp(8)+comp(9)+comp(10))/nbz, !O/Z 3 (comp(11)+comp(12)+comp(13))/nbz, !Ne/Z 4 comp(14)/nbz, !Na/Z 5 (comp(15)+comp(16)+comp(17))/nbz, !Mg/Z 6 comp(18)/nbz, !Al/Z 7 (comp(19)+comp(20))/nbz, !Si/Z 8 comp(21)/nbz, !n/Z 9 comp(22)/nbz !Ex/Z WRITE(2,6) ; WRITE(*,6) 6 FORMAT(/,'H2, Li7, Be7 à l''équilibre') 21 FORMAT(/,'Rapports en nombre dans Z:',/,'C/Z:',es10.3,', N/Z:',es10.3, 1 ', O/Z:',es10.3,', Ne/Z:',es10.3,', Na/Z:',es10.3,/, 2 'Mg/Z:',es10.3,', Al/Z:',es10.3,', Si/Z:',es10.3, 3 ', n/Z:',es10.3,/,'Ex/Z:',es10.3) WRITE(2,14)he3she4z,c13sc12,n15sn14,o17so16,o18so16,ne21sne20,ne22sne20 WRITE(*,14)he3she4z,c13sc12,n15sn14,o17so16,o18so16,ne21sne20,ne22sne20 14 FORMAT(/,'Rapports isotopiques en nombre:',/, 1 'He3/He4=',es10.3,', C13/C12=',es10.3, 2 ', N15/N14=',es10.3,/,'O17/O16=',es10.3,', O18/O16=',es10.3, 3',Ne21/Ne20=',es10.3,',Ne22/Ne20=',es10.3) WRITE(*,15)mg25smg24,mg26smg25,si29ssi28 WRITE(2,15)mg25smg24,mg26smg25,si29ssi28 15 FORMAT('Mg25/Mg24=',es10.3,', Mg26Mg25=',es10.3,', Si29/Si28=',es10.3) WRITE(2,7) ; WRITE(*,7) 7 FORMAT(/,'on utilise une table') c définitions diverses anuc=ANINT(nucleo) !nombre atomique c nettoyage DEALLOCATE(a,b,comp_dex,indpc) c------------------------------------------------------------------ c les réactions CASE(2) dcomp=0.d0 ; jac=0.d0 IF(t < t_inf)RETURN c WRITE(*,2000)t,ro ; PRINT*,'entrée' CALL rq_reac(comp,t,ro,r,drt,dro,drx,q,dqt,dqo,dqx,mue,dmuex) c WRITE(*,*)'comp' ; WRITE(*,2000)comp(1:nchim) c WRITE(*,*)'réactions' ; WRITE(*,2000)r(1:nreac) c test de suppression des neutrons c r(22:33)=0.d0 ; drx(22:33,:)=0.d0 c équations d'évolution c équation dcomp(H1) 1 dcomp(1)=-(3.d0*r(1)*comp(1)+r(8)*comp(4)+r(9)*comp(5)+r(10)*comp(6) 1 +(r(11)+r(12))*comp(7)+r(13)*comp(8)+r(14)*comp(9)+r(34)*comp(10) 2 +r(35)*comp(9)+(r(36)+r(37))*comp(18)+r(38)*comp(14) 3 +r(39)*comp(13))*comp(1)+(2.d0*r(3)*comp(2)-r(4)*comp(3))*comp(2) 4 +r(45)*comp(4)**2+r(48)*comp(8)*comp(4) c équation dcomp(He3) 2 dcomp(2)=r(1)*comp(1)**2-(2.d0*r(3)*comp(2)+r(4)*comp(3))*comp(2) c équation dcomp(He4) 3 dcomp(3)=(r(3)*comp(2)+r(4)*comp(3))*comp(2) 1 +(r(12)*comp(7)+r(14)*comp(9)+r(34)*comp(10)+r(36)*comp(18))*comp(1) 2 -(3.d0*r(15)*comp(3)**2+r(16)*comp(4)+r(17)*comp(8)+r(18)*comp(6) 3 +r(19)*comp(10)+(r(20)+r(21))*comp(13)+r(25)*comp(12)+r(31)*comp(5) 4 +r(32)*comp(16)+r(33)*comp(17)+r(42)*comp(11)+r(43)*comp(15))*comp(3) 5 +r(40)*comp(11)+r(41)*comp(15)+r(46)*comp(4)**2+r(49)*comp(8)*comp(4) 6 +r(50)*comp(8)**2 c équation dcomp(C12) 4 dcomp(4)=(-r(8)*comp(4)+r(12)*comp(7))*comp(1) 1 +(r(15)*comp(3)**2-r(16)*comp(4))*comp(3) 2 -(r(44)+r(45)+r(46)+r(51))*comp(4)**2-(r(47)+r(48)+r(49))*comp(8)*comp(4) c équation dcomp(C13) 5 dcomp(5)=(r(8)*comp(4)-r(9)*comp(5))*comp(1)-r(31)*comp(5)*comp(3) c équation dcomp(N14) 6 dcomp(6)=(r(9)*comp(5)-r(10)*comp(6)+r(14)*comp(9))*comp(1) 1 -r(18)*comp(6)*comp(3) c équation dcomp(N15) 7 dcomp(7)=(r(10)*comp(6)-(r(11)+r(12))*comp(7)+r(34)*comp(10))*comp(1) c équation dcomp(O16) 8 dcomp(8)=(r(11)*comp(7)-r(13)*comp(8))*comp(1) 1 +(r(16)*comp(4)-r(17)*comp(8)+r(31)*comp(5))*comp(3)+r(40)*comp(11) 2 -(r(47)+r(48)+r(49))*comp(8)*comp(4)-r(50)*comp(8)**2 c équation dcomp(O17) 9 dcomp(9)=(r(13)*comp(8)-r(14)*comp(9)-r(35)*comp(9))*comp(1) c équation dcomp(O18) 10 dcomp(10)=r(18)*comp(6)*comp(3)-r(19)*comp(10)*comp(3) 1 -r(34)*comp(10)*comp(1)+r(35)*comp(9)*comp(1) c équation dcomp(Ne20) 11 dcomp(11)=r(17)*comp(8)*comp(3)-r(24)*comp(11)*comp(21)-r(40)*comp(11) 1 +r(41)*comp(15)-r(42)*comp(11)*comp(3)+r(46)*comp(4)**2 c équation dcomp(Ne21) 12 dcomp(12)=r(24)*comp(11)*comp(21)-r(25)*comp(12)*comp(3) c équation dcomp(Ne22) 13 dcomp(13)=r(19)*comp(10)*comp(3)-(r(20)+r(21))*comp(13)*comp(3) 1 -r(22)*comp(13)*comp(21)-r(39)*comp(13)*comp(1) c équation dcomp(Na23) 14 dcomp(14)=r(22)*comp(13)*comp(21)-r(23)*comp(14)*comp(21) 1 -r(38)*comp(14)*comp(1)+r(39)*comp(13)*comp(1) 2 +(r(45)+r(51))*comp(4)**2 c équation dcomp(Mg24) 15 dcomp(15)=r(23)*comp(14)*comp(21)+r(25)*comp(12)*comp(3) 1 -r(26)*comp(15)*comp(21)+(r(36)*comp(18)+r(38)*comp(14))*comp(1) 2 -r(41)*comp(15)+r(42)*comp(11)*comp(3)-r(43)*comp(15)*comp(3) 3 +r(44)*comp(4)**2+r(49)*comp(8)*comp(4) c équation dcomp(Mg25) 16 dcomp(16)=r(21)*comp(13)*comp(3)+r(26)*comp(15)*comp(21) 1 -r(27)*comp(16)*comp(21)-r(32)*comp(16)*comp(3) c équation dcomp(Mg26) 17 dcomp(17)=r(20)*comp(13)*comp(3)+r(27)*comp(16)*comp(21) 1 -r(28)*comp(17)*comp(21)-r(33)*comp(17)*comp(3) c équation dcomp(Al27) 18 dcomp(18)=r(28)*comp(17)*comp(21)-r(29)*comp(18)*comp(21) 1 -(r(36)+r(37))*comp(18)*comp(1)+r(48)*comp(8)*comp(4) c équation dcomp(Si28) 19 dcomp(19)=r(29)*comp(18)*comp(21)-r(30)*comp(19)*comp(21) 1 +r(32)*comp(16)*comp(3)+r(37)*comp(18)*comp(1)+r(43)*comp(15)*comp(3) 2 +r(47)*comp(8)*comp(4)+r(50)*comp(8)**2 c équation dcomp(Si29) 20 dcomp(20)=r(30)*comp(19)*comp(21)+r(33)*comp(17)*comp(3) c équation dcomp(n) 21 dcomp(21)=(r(21)*comp(13)+r(25)*comp(12)+r(31)*comp(5) 1 +r(32)*comp(16)+r(33)*comp(17))*comp(3) 2 -(r(22)*comp(13)+r(23)*comp(14)+r(24)*comp(11)+r(26)*comp(15) 3 +r(27)*comp(16)+r(28)*comp(17)+r(29)*comp(18)+r(30)*comp(19))*comp(21) 4 +r(51)*comp(4)**2 c équation dcomp(Ex=0) i_ex=22 tous les isotopes étant pris en compte, c les équations assurent la conservation des nucléons c Sum(X_i)=1 est assuré dans evol en normalisant dcomp(i_ex)=-DOT_PRODUCT(dcomp,anuc)/anuc(i_ex) c calcul du jacobien IF(deriv)THEN !jac(i,j) : équation, j : élément i c équation dcomp(1), H1 c dcomp(1)=-(3.d0*r(1)*comp(1)+r(8)*comp(4)+r(9)*comp(5)+r(10)*comp(6)+ c 1 (r(11)+r(12))*comp(7)+r(13)*comp(8)+r(14)*comp(9)+r(34)*comp(10) c 2 +r(35)*comp(9)+(r(36)+r(37))*comp(18)+r(38)*comp(14) c 3 +r(39)*comp(13))*comp(1)+(2.d0*r(3)*comp(2)-r(4)*comp(3))*comp(2) c 4 +r(45)*comp(4)**2+r(48)*comp(8)*comp(4) jac(1,1)=-6.d0*r(1)*comp(1)-r(8)*comp(4)-r(9)*comp(5)-r(10)*comp(6) 1 -(r(11)+r(12))*comp(7)-r(13)*comp(8)-r(14)*comp(9) 2 -r(34)*comp(10)-r(35)*comp(9)-(r(36)+r(37))*comp(18) 3 -r(38)*comp(14)-r(39)*comp(13) !d /H1 (1) jac(1,2)=4.d0*r(3)*comp(2)-r(4)*comp(3) !d /He3 (2) jac(1,3)=-r(4)*comp(2) !d /He4 (3) jac(1,4)=-r(8)*comp(1)+2.d0*r(45)*comp(4) 1 +r(48)*comp(8) !d /C12 (4) jac(1,5)=-r(9)*comp(1) !d /C13 (5) jac(1,6)=-r(10)*comp(1) !d /N14 (6) jac(1,7)=-(r(11)+r(12))*comp(1) !d /N15 (7) jac(1,8)=-r(13)*comp(1)+r(48)*comp(4) !d /O16 (8) jac(1,9)=-(r(14)+r(35))*comp(1) !d /O17 (9) jac(1,10)=-r(34)*comp(1) !d /O18 (10) jac(1,13)=-r(39)*comp(1) !d /Ne22 (13) jac(1,14)=-r(38)*comp(1) !d /Na23 (14) jac(1,18)=-(r(36)+r(37))*comp(1) !d /Al27 (18) DO i=1,nchim !dépendances dues à l'effet d'écran jac(1,i)=jac(1,i) 1 -(3.d0*drx(1,i)*comp(1)+drx(8,i)*comp(4)+drx(9,i)*comp(5) 2 +drx(10,i)*comp(6)+(drx(11,i)+drx(12,i))*comp(7)+drx(13,i)*comp(8) 2 +drx(14,i)*comp(9)+drx(34,i)*comp(10)+drx(35,i)*comp(9) 3 +(drx(36,i)+drx(37,i))*comp(18)+drx(38,i)*comp(14))*comp(1) 4 +(2.d0*drx(3,i)*comp(2)-drx(4,i)*comp(3))*comp(2)+drx(45,i)*comp(4)**2 5 +drx(48,i)*comp(8)*comp(4) ENDDO c équation dcomp(He3) 2 c dcomp(2)=r(1)*comp(1)**2-(2.d0*r(3)*comp(2)+r(4)*comp(3))*comp(2) jac(2,1)=2.d0*r(1)*comp(1) !d /H1 (1) jac(2,2)=-4.d0*r(3)*comp(2)-r(4)*comp(3) !d /He3 (2) jac(2,3)=-r(4)*comp(2) !d /He4 (3) DO i=1,nchim !dépendances dues à l'effet d'écran jac(2,i)=jac(2,i) 1 +drx(1,i)*comp(1)**2-(2.d0*drx(3,i)*comp(2)+drx(4,i)*comp(3))*comp(2) ENDDO c équation dcomp(He4) 3 c dcomp(3)=(r(3)*comp(2)+r(4)*comp(3))*comp(2) c 1 +(r(12)*comp(7)+r(14)*comp(9)+r(34)*comp(10)+r(36)*comp(18))*comp(1) c 2 -(3.d0*r(15)*comp(3)**2+r(16)*comp(4)+r(17)*comp(8)+r(18)*comp(6) c 3 +r(19)*comp(10)+(r(20)+r(21))*comp(13)+r(25)*comp(12)+r(31)*comp(5) c 4 +r(32)*comp(16)+r(33)*comp(17))*comp(3)+r(40)*comp(11)+r(41)*comp(15) c 5 -r(42)*comp(11)*comp(3)-r(43)*comp(15)*comp(3)+r(46)*comp(4)**2 c 6 +r(49)*comp(8)*comp(4)+r(50)*comp(8)**2 jac(3,1)=r(12)*comp(7)+r(14)*comp(9)+r(34)*comp(10) 1 +r(36)*comp(18) !d /H1 (1) jac(3,2)=2.d0*r(3)*comp(2)+r(4)*comp(3) !d /He3 (2) jac(3,3)=r(4)*comp(2)-9.d0*r(15)*comp(3)**2-r(16)*comp(4) 1 -r(17)*comp(8)-r(18)*comp(6)-r(19)*comp(10)-(r(20)+r(21))*comp(13) 2 -r(25)*comp(12)-r(31)*comp(5)-r(32)*comp(16)-r(33)*comp(17) 3 -r(42)*comp(11)-r(43)*comp(15) !d /He4 (3) jac(3,4)=-r(16)*comp(3)+2.d0*r(46)*comp(4)+r(49)*comp(8) !d /C12 (4) jac(3,5)=-r(31)*comp(3) !d /C13 (5) jac(3,6)=-r(18)*comp(3) !d /N14 (6) jac(3,7)=r(12)*comp(1) !d /N15 (7) jac(3,8)=-r(17)*comp(3)+r(49)*comp(4)+2.d0*r(50)*comp(8) !d /O16 (8) jac(3,9)=r(14)*comp(1) !d /O17 (9) jac(3,10)=-r(19)*comp(3)+r(34)*comp(1) !d /O18 (10) jac(3,11)=r(40)-r(42)*comp(3) !d /Ne20 (11) jac(3,12)=-r(25)*comp(3) !d /Ne21 (12) jac(3,13)=-(r(20)+r(21))*comp(3) !d /Ne22 (13) jac(3,15)=r(41)-r(43)*comp(3) !d /Mg24 (15) jac(3,16)=-r(32)*comp(3) !d /Mg25 (16) jac(3,17)=-r(33)*comp(3) !d /Mg26 (17) jac(3,18)=r(36)*comp(1) !d /Al27 (18) DO i=1,nchim !dépendances dues à l'effet d'écran jac(3,i)=jac(3,i) 1 +(drx(3,i)*comp(2)+drx(4,i)*comp(3))*comp(2) 2 +(drx(12,i)*comp(7)+drx(14,i)*comp(9)+drx(34,i)*comp(10))*comp(1) 3 -(3.d0*drx(15,i)*comp(3)**2+drx(16,i)*comp(4)+drx(17,i)*comp(8) 4 +drx(18,i)*comp(6)+drx(19,i)*comp(10)+(drx(20,i)+drx(21,i))*comp(13) 5 +drx(25,i)*comp(12)+drx(31,i)*comp(5)+drx(32,i)*comp(16) 6 +drx(33,i)*comp(17))*comp(3)+drx(40,i)*comp(11)+drx(41,i)*comp(15) 7 -drx(42,i)*comp(11)*comp(3)-drx(43,i)*comp(15)*comp(3) 8 +drx(46,i)*comp(4)**2+drx(49,i)*comp(8)*comp(4)+drx(50,i)*comp(8)**2 ENDDO c équation dcomp(C12) 4 c dcomp(4)=(-r(8)*comp(4)+r(12)*comp(7))*comp(1) c 1 +(r(15)*comp(3)**2-r(16)*comp(4))*comp(3) c 2 -(r(44)+r(45)+r(46)+r(51))*comp(4)**2-(r(47)+r(48)+r(49))*comp(8)*comp(4) jac(4,1)=-r(8)*comp(4)+r(12)*comp(7) !d /H1 (1) jac(4,3)=3.d0*r(15)*comp(3)**2-r(16)*comp(4) !d /He4 (3) jac(4,4)=-r(8)*comp(1)-r(16)*comp(3) 1 -2.d0*(r(44)+r(45)+r(46)+r(51))*comp(4) 2 -(r(47)+r(48)+r(49))*comp(8) !d /C12 (4) jac(4,7)=r(12)*comp(1) !d /N15 (7) jac(4,8)=-(r(47)+r(48)+r(49))*comp(4) !d /O16 (8) DO i=1,nchim !dépendances dues à l'effet d'écran jac(4,i)=jac(4,i) 1 +(-drx(8,i)*comp(4)+drx(12,i)*comp(7))*comp(1) 2 +(drx(15,i)*comp(3)**2-drx(16,i)*comp(4))*comp(3) 3 -(drx(44,i)+drx(45,i)+drx(46,i)+drx(51,i))*comp(4)**2 4 -(drx(47,i)+drx(48,i)+drx(49,i))*comp(8)*comp(4) ENDDO c équation dcomp(C13) 5 c dcomp(5)=(r(8)*comp(4)-r(9)*comp(5))*comp(1)-r(31)*comp(5)*comp(3) jac(5,1)=r(8)*comp(4)-r(9)*comp(5) !d /H1 (1) jac(5,3)=-r(31)*comp(5) !d /He4 (3) jac(5,4)=r(8)*comp(1) !d /C12 (4) jac(5,5)=-r(9)*comp(1)-r(31)*comp(3) !d /C13 (5) DO i=1,nchim !dépendances dues à l'effet d'écran jac(5,i)=jac(5,i) 1 +(drx(8,i)*comp(4)-drx(9,i)*comp(5))*comp(1)-drx(31,i)*comp(5)*comp(3) ENDDO c équation dcomp(N14) 6 c dcomp(6)=(r(9)*comp(5)-r(10)*comp(6)+r(14)*comp(9))*comp(1) c -r(18)*comp(6)*comp(3) jac(6,1)=r(9)*comp(5)-r(10)*comp(6)+r(14)*comp(9) !d /H1 (1) jac(6,3)=-r(18)*comp(6) !d /He4 (3) jac(6,5)=r(9)*comp(1) !d /C13 (5) jac(6,6)=-r(10)*comp(1)-r(18)*comp(3) !d /N14 (6) jac(6,9)=r(14)*comp(1) !d /O17 (9) DO i=1,nchim !dépendances dues à l'effet d'écran jac(6,i)=jac(6,i) 1 +(drx(9,i)*comp(5)-drx(10,i)*comp(6)+drx(14,i)*comp(9))*comp(1) 2 -drx(18,i)*comp(6)*comp(3) ENDDO c équation dcomp(N15) 7 c dcomp(7)=(r(10)*comp(6)-(r(11)+r(12))*comp(7)+r(34)*comp(10))*comp(1) jac(7,1)=r(10)*comp(6)-(r(11)+r(12))*comp(7)+r(34)*comp(10) !d /H1 (1) jac(7,6)=r(10)*comp(1) !d /N14 (6) jac(7,7)=-(r(11)+r(12))*comp(1) !d /N15 (7) jac(7,10)=r(34)*comp(1) !d /O18 (10) DO i=1,nchim !dépendances dues à l'effet d'écran jac(7,i)=jac(7,i) 1 +(drx(10,i)*comp(6)-(drx(11,i)+drx(12,i))*comp(7) 2 +drx(34,i)*comp(10))*comp(1) ENDDO c équation dcomp(O16) 8 c dcomp(8)=(r(11)*comp(7)-r(13)*comp(8))*comp(1) c 1 +(r(16)*comp(4)-r(17)*comp(8)+r(31)*comp(5))*comp(3)+r(40)*comp(11) c 2 -(r(47)+r(48)+r(49))*comp(8)*comp(4)-r(50)*comp(8)**2 jac(8,1)=r(11)*comp(7)-r(13)*comp(8) !d /H1 (1) jac(8,3)=r(16)*comp(4)-r(17)*comp(8)+r(31)*comp(5) !d /He4 (3) jac(8,4)=r(16)*comp(3)-(r(47)+r(48)+r(49))*comp(8) !d /C12 (4) jac(8,5)=r(31)*comp(3) !d /C13 (5) jac(8,7)=r(11)*comp(1) !d /N15 (7) jac(8,8)=-r(13)*comp(1)-r(17)*comp(3)-(r(47)+r(48)+r(49))*comp(4) 1 -2.d0*r(50)*comp(8) !d /O16 (8) jac(8,11)=r(40) !d /Ne20 (11) DO i=1,nchim !dépendances dues à l'effet d'écran jac(8,i)=jac(8,i) 1 +(drx(11,i)*comp(7)-drx(13,i)*comp(8))*comp(1) 2 +(drx(16,i)*comp(4)-drx(17,i)*comp(8)+drx(31,i)*comp(5))*comp(3) 3 +drx(40,i)*comp(11)-(drx(47,i)+drx(48,i)+drx(49,i))*comp(8)*comp(4) 4 -drx(50,i)*comp(8)**2 ENDDO c équation dcomp(O17) 9 c dcomp(9)=(r(13)*comp(8)-r(14)*comp(9)-r(35)*comp(9))*comp(1) jac(9,1)=r(13)*comp(8)-r(14)*comp(9)-r(35)*comp(9) !d /H1 (1) jac(9,8)=r(13)*comp(1) !d /O16 (8) jac(9,9)=-r(14)*comp(1)-r(35)*comp(1) !d /O17 (9) DO i=1,nchim !dépendances dues à l'effet d'écran jac(9,i)=jac(9,i) 1 +(drx(13,i)*comp(8)-drx(14,i)*comp(9)-drx(35,i)*comp(9))*comp(1) ENDDO c équation dcomp(O18) 10 c dcomp(10)=r(18)*comp(6)*comp(3)-r(19)*comp(10)*comp(3) c 1 -r(34)*comp(10)*comp(1)+r(35)*comp(9)*comp(1) jac(10,1)=-r(34)*comp(10)+r(35)*comp(9) !d /H1 (1) jac(10,3)=r(18)*comp(6)-r(19)*comp(10) !d /He4 (3) jac(10,6)=r(18)*comp(3) !d /N14 (6) jac(10,9)=r(35)*comp(1) !d /O17 (9) jac(10,10)=-r(19)*comp(3)-r(34)*comp(1) !d /O18 (10) DO i=1,nchim !dépendances dues à l'effet d'écran jac(10,i)=jac(10,i) 1 +drx(18,i)*comp(6)*comp(3)-drx(19,i)*comp(10)*comp(3) 2 -drx(34,i)*comp(10)*comp(1)+drx(35,i)*comp(9)*comp(1) ENDDO c équation dcomp(Ne20) 11 c dcomp(11)=r(17)*comp(8)*comp(3)-r(24)*comp(11)*comp(21) c 1 -r(40)*comp(11)+r(41)*comp(15)-r(42)*comp(11)*comp(3) c 2 -r(42)*comp(11)*comp(3)+r(46)*comp(4)**2 jac(11,3)=r(17)*comp(8)-r(42)*comp(11) !d /He4 (3) jac(11,8)=r(17)*comp(3) !d /O16 (8) jac(11,11)=-r(24)*comp(21)-r(40)-r(42)*comp(3) 1 +2.d0*r(46)*comp(4) !d /N20 (11) jac(11,15)=r(41) !d /Mg24 (15) jac(11,21)=-r(24)*comp(11) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(11,i)=jac(11,i) 1 +drx(17,i)*comp(8)*comp(3)-drx(24,i)*comp(11)*comp(21) 2 -drx(40,i)*comp(11)+drx(41,i)*comp(15)-drx(42,i)*comp(11)*comp(3) 3 +drx(46,i)*comp(4)**2 ENDDO c équation dcomp(Ne21) 12 c dcomp(12)=r(24)*comp(11)*comp(21)-r(25)*comp(12)*comp(3) jac(12,3)=-r(25)*comp(12) !d /He4 (3) jac(12,11)=r(24)*comp(21) !d /Ne20 (11) jac(12,12)=-r(25)*comp(3) !d /Ne21 (12) jac(12,21)=r(24)*comp(11) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(12,i)=jac(12,i) 1 +drx(24,i)*comp(11)*comp(21)-drx(25,i)*comp(12)*comp(3) ENDDO c équation dcomp(Ne22) 13 c dcomp(13)=r(19)*comp(10)*comp(3)-(r(20)+r(21))*comp(13)*comp(3) c 1 -r(22)*comp(13)*comp(21)-r(39)*comp(13)*comp(1) jac(13,1)=-r(39)*comp(13) !d /H1 (1) jac(13,3)=r(19)*comp(10)-(r(20)+r(21))*comp(13) !d /He4 (3) jac(13,10)=r(19)*comp(3) !d /O18 (10) jac(13,13)=-(r(20)+r(21))*comp(3)-r(39)*comp(1) !d /Ne22 (13) jac(13,21)=-r(22)*comp(13) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(13,i)=jac(13,i) 1 +drx(19,i)*comp(10)*comp(3)-(drx(20,i)+drx(21,i))*comp(13)*comp(3) 2 -drx(22,i)*comp(13)*comp(21)-drx(39,i)*comp(13)*comp(1) ENDDO c équation dcomp(Na23) 14 c dcomp(14)=r(22)*comp(13)*comp(21)-r(23)*comp(14)*comp(21) c 1 -r(38)*comp(14)*comp(1)+r(39)*comp(13)*comp(1)+(r(45)+r(51))*comp(4)**2 jac(14,1)=-r(38)*comp(14)+r(39)*comp(13) !d /H1 (1) jac(14,4)=2.d0*(r(45)+r(51))*comp(4) !d /C12 (4) jac(14,13)=r(22)*comp(21)+r(39)*comp(1) !d /Ne22 (13) jac(14,14)=-r(23)*comp(21)-r(38)*comp(1) !d /Na23 (14) jac(14,21)=r(22)*comp(13)-r(23)*comp(14) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(14,i)=jac(14,i) 1 +drx(22,i)*comp(13)*comp(21)-drx(23,i)*comp(14)*comp(21) 2 +(-drx(38,i)*comp(14)+drx(39,i)*comp(13))*comp(1) 3 +(drx(45,i)+drx(51,i))*comp(4)**2 ENDDO c équation dcomp(Mg24) 15 c dcomp(15)=r(23)*comp(14)*comp(21)+r(25)*comp(12)*comp(3) c 1 -r(26)*comp(15)*comp(21)+r(36)*comp(18)*comp(1) c 2 +r(38)*comp(14)*comp(1)-r(41)*comp(15)+r(42)*comp(11)*comp(3) c 3 +r(42)*comp(11)*comp(3)-r(43)*comp(15)*comp(3)+r(44)*comp(4)**2 c 4 +r(49)*comp(8)*comp(4) jac(15,1)=r(36)*comp(18)+r(38)*comp(14) !d /H1 (1) jac(15,3)=r(25)*comp(12)+r(42)*comp(11)-r(43)*comp(15)!d /He4 (3) jac(15,4)=2.d0*r(44)*comp(4)+r(49)*comp(8) !d /C12 (4) jac(15,8)=r(49)*comp(4) !d /O16 jac(15,11)=r(42)*comp(3) !d /Ne20 (11) jac(15,12)=r(25)*comp(3) !d /Ne21 (12) jac(15,14)=r(23)*comp(21)+r(38)*comp(1) !d /Na23 (14) jac(15,15)=-r(26)*comp(21)-r(41)-r(43)*comp(3) !d /Mg24 (15) jac(15,18)=r(36)*comp(1) !d /Al27 (18) jac(15,21)=r(23)*comp(14)-r(26)*comp(15) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(15,i)=jac(15,i) 1 +drx(23,i)*comp(14)*comp(21)+drx(25,i)*comp(12)*comp(3) 2 -drx(26,i)*comp(15)*comp(21)+(drx(36,i)*comp(18) 3 +drx(38,i)*comp(14))*comp(1)-drx(41,i)*comp(15) 4 +drx(42,i)*comp(11)*comp(3)-drx(43,i)*comp(15)*comp(3) 5 +drx(44,i)*comp(4)**2+drx(49,i)*comp(8)*comp(4) ENDDO c équation dcomp(Mg25) 16 c dcomp(16)=r(21)*comp(13)*comp(3)+r(26)*comp(15)*comp(21) c 1 -r(27)*comp(16)*comp(21)-r(32)*comp(16)*comp(3) jac(16,3)=r(21)*comp(13)-r(32)*comp(16) !d /He4 (3) jac(16,13)=r(21)*comp(3) !d /Ne22 (13) jac(16,15)=r(26)*comp(21) !d /Mg24 (15) jac(16,16)=-r(27)*comp(21)-r(32)*comp(3) !d /Mg25 (16) jac(16,21)=r(26)*comp(15)-r(27)*comp(16) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(16,i)=jac(16,i) 1 +drx(21,i)*comp(13)*comp(3)+drx(26,i)*comp(15)*comp(21) 2 -drx(27,i)*comp(16)*comp(21)-drx(32,i)*comp(16)*comp(3) ENDDO c équation dcomp(Mg26) 17 c dcomp(17)=r(20)*comp(13)*comp(3)+r(27)*comp(16)*comp(21) c 1 -r(28)*comp(17)*comp(21)-r(33)*comp(17)*comp(3) jac(17,3)=r(20)*comp(13)-r(33)*comp(17) !d /He4 (3) jac(17,13)=r(20)*comp(3) !d /Ne22 (13) jac(17,16)=r(27)*comp(21) !d /Mg25 (16) jac(17,17)=-r(28)*comp(21)-r(33)*comp(3) !d /Mg26 (17) jac(17,21)=r(27)*comp(16)-r(28)*comp(17) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(17,i)=jac(17,i) 1 +drx(20,i)*comp(13)*comp(3)+drx(27,i)*comp(16)*comp(21) 2 -drx(28,i)*comp(17)*comp(21)-drx(33,i)*comp(17)*comp(3) ENDDO c équation dcomp(Al27) 18 c dcomp(18)=r(28)*comp(17)*comp(21)-r(29)*comp(18)*comp(21) c 1 -(r(36)+r(37))*comp(18)*comp(1)+r(48)*comp(8)*comp(4) jac(18,1)=-(r(36)+r(37))*comp(18) !d /H1 (1) jac(18,4)=r(48)*comp(8) !d /C12 jac(18,8)=r(48)*comp(4) !d /O16 jac(18,17)=r(28)*comp(21) !d /Mg26 (17) jac(18,18)=-r(29)*comp(21)-(r(36)+r(37))*comp(1) !d /Al27 (18) jac(18,21)=r(28)*comp(17)-r(29)*comp(18) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(18,i)=jac(18,i) 1 +drx(28,i)*comp(17)*comp(21)-drx(29,i)*comp(18)*comp(21) 2 -(drx(36,i)+drx(37,i))*comp(18)*comp(1)+drx(48,i)*comp(8)*comp(4) ENDDO c équation dcomp(Si28) 19 c dcomp(19)=r(29)*comp(18)*comp(21)-r(30)*comp(19)*comp(21) c 1 +r(32)*comp(16)*comp(3)+r(37)*comp(18)*comp(1)+r(43)*comp(15)*comp(3) c 2 +r(47)*comp(8)*comp(4)+r(50)*comp(8)**2 jac(19,1)=r(37)*comp(18) !d /H1 (1) jac(19,3)=r(32)*comp(16)+r(43)*comp(15)+r(47)*comp(8) !d /He4 (3) jac(19,8)=r(47)*comp(4)+2.d0*r(50)*comp(8) !d /O16 (8) jac(19,15)=r(43)*comp(3) !d /Mg24 (15) jac(19,16)=r(32)*comp(3) !d /Mg25 (16) jac(19,18)=r(29)*comp(21)+r(37)*comp(1) !d /Al27 (18) jac(19,19)=-r(30)*comp(21) !d /Si28 (19) jac(19,21)=r(29)*comp(18)-r(30)*comp(19) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(19,i)=jac(19,i) 1 +drx(29,i)*comp(18)*comp(21)-drx(30,i)*comp(19)*comp(21) 2 +drx(32,i)*comp(16)*comp(3)+drx(37,i)*comp(18)*comp(1) 3 +drx(43,i)*comp(15)*comp(3)+drx(47,i)*comp(8)*comp(4) 4 +drx(50,i)*comp(8)**2 ENDDO c équation dcomp(Si29) 20 c dcomp(20)=r(30)*comp(19)*comp(21)+r(33)*comp(17)*comp(3) jac(20,3)=r(33)*comp(17) !d /He4 (3) jac(20,17)=r(33)*comp(3) !d /Mg26 (17) jac(20,19)=r(30)*comp(21) !d /Si28 (19) jac(20,21)=r(30)*comp(19) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(20,i)=jac(20,i) 1 +drx(30,i)*comp(19)*comp(21)+drx(33,i)*comp(17)*comp(3) ENDDO c équation dcomp(n) 21 c dcomp(21)=(r(21)*comp(13)+r(25)*comp(12)+r(31)*comp(5) c 1 +r(32)*comp(16)+r(33)*comp(17))*comp(3) c 2 -(r(22)*comp(13)+r(23)*comp(14)+r(24)*comp(11)+r(26)*comp(15) c 3 +r(27)*comp(16)+r(28)*comp(17)+r(29)*comp(18)+r(30)*comp(19))*comp(21) c 4 +r(51)*comp(4)**2 jac(21,3)=r(21)*comp(13)+r(25)*comp(12)+r(31)*comp(5) 1 +r(32)*comp(16)+r(33)*comp(17) !d /He4 (3) jac(21,4)=2.d0*r(51)*comp(4) !d /C12 (4) jac(21,5)=r(31)*comp(3) !d /C13 (5) jac(21,11)=-r(24)*comp(21) !d /Ne20 (11) jac(21,12)=r(25)*comp(3) !d /Ne21 (12) jac(21,13)=r(21)*comp(3)-r(22)*comp(21) !d /Ne22 (13) jac(21,14)=-r(23)*comp(21) !d /Na23 (14) jac(21,15)=-r(26)*comp(21) !d /Mg24 (15) jac(21,16)=r(32)*comp(3)-r(27)*comp(21) !d /Mg25 (16) jac(21,17)=r(33)*comp(3)-r(28)*comp(21) !d /Mg26 (17) jac(21,18)=-r(29)*comp(21) !d /Al27 (18) jac(21,19)=-r(30)*comp(21) !d /Si28 (19) jac(21,21)=-r(22)*comp(13)-r(23)*comp(14)-r(24)*comp(11) 1 -r(26)*comp(15)-r(27)*comp(16)-r(28)*comp(17)+r(25)*comp(12) 2 -r(29)*comp(18)-r(30)*comp(19) !d /n (21) DO i=1,nchim !dépendances dues à l'effet d'écran jac(21,i)=jac(21,i) 1 +(drx(21,i)*comp(13)+drx(25,i)*comp(12)+drx(31,i)*comp(5) 2 +drx(32,i)*comp(16)+drx(33,i)*comp(17))*comp(3) 3 -(drx(22,i)*comp(13)+drx(23,i)*comp(14)+drx(24,i)*comp(11) 4 +drx(26,i)*comp(15)+drx(27,i)*comp(16) 5 +drx(28,i)*comp(17)+drx(29,i)*comp(18)+drx(30,i)*comp(19))*comp(21) 6 +drx(51,i)*comp(4)**2 ENDDO c équation dcomp(Ex=0) i_ex=22 tous les isotopes étant pris en compte, c les équations assurent la conservation des nucléons c Sum(X_i)=1 est assuré dans evol en normalisant c dcomp(i_ex)=-DOT_PRODUCT(dcomp,anuc)/anuc(i_ex) DO j=1,nchim B2: DO i=1,nchim IF(i == i_ex)CYCLE B2 jac(i_ex,j)=jac(i_ex,j)+anuc(i)*jac(i,j) ENDDO B2 ENDDO jac(i_ex,:)=-jac(i_ex,:)/anuc(i_ex) c unités de temps pour intégration temporelle jac=jac*secon6 ENDIF !deriv c unités de temps pour intégration temporelle dcomp=dcomp*secon6 c------------------------------------------------------------------- CASE(3) c calcul de la production d'énergie nucléaire et dérivées c pour H2(H,g)He3, q(2)H**2=q(2)*r(1)/r(2) epsilon(1:4)=0.d0 ; et=0.d0 ; ero=0.d0 ; ex=0.d0 IF(t <= t_inf)RETURN CALL rq_reac(comp,t,ro,r,drt,dro,drx,q,dqt,dqo,dqx,mue,dmuex) c PRINT*,'les q' ; WRITE(*,2000)q c test de suppression des neutrons c q(22:33)=0.d0 ; dqt(22:33)=0.d0 ; dqo(22:33)=0.d0 ; dqx(22:33,:)=0.d0 c mue : nombre d'électrons / mole /g = 1/poids mol. moy. par e- IF(comp(1) > 0.d0)THEN h2=r(1)/r(2)*comp(1) ; den=r(6)*mue+r(7)*comp(1) be7=r(4)*comp(2)*comp(3)/den ; li7=r(6)*be7*mue/r(5)/comp(1) ELSE h2=0.d0 ; be7=0.d0 ; li7=0.d0 ENDIF c PRINT*,'h2,li7,be7' ; WRITE(*,2000)h2,li7,be7 c énergie PP epsilon(2)=(q(1)*comp(1)+q(2)*h2+q(5)*li7+q(7)*be7)*comp(1) 1 +(q(3)*comp(2)+q(4)*comp(3))*comp(2)+q(6)*mue*be7 c énergie CNO epsilon(3)=(q(8)*comp(4)+q(9)*comp(5)+q(10)*comp(6) 1 +(q(11)+q(12))*comp(7)+q(13)*comp(8)+q(14)*comp(9)+q(34)*comp(10) 2 +q(35)*comp(9)+(q(36)+q(37))*comp(18)+q(38)*comp(14) 3 +r(39)*comp(13))*comp(1) c énergies 3alpha C Ne O epsilon(4)=(q(15)*comp(3)**2+q(16)*comp(4)+q(17)*comp(8) 1 +q(18)*comp(6)+q(19)*comp(10)+(q(20)+q(21))*comp(13)+q(25)*comp(12) 2 +q(31)*comp(5)+q(32)*comp(16)+q(33)*comp(17)+q(42)*comp(11) 3 +q(43)*comp(15))*comp(3)+(q(44)+q(45)+q(46)+q(51))*comp(4)**2 4 +(q(22)*comp(13)+q(23)*comp(14)+q(24)*comp(11)+q(26)*comp(15) 5 +q(27)*comp(16)+q(28)*comp(17)+q(29)*comp(18)+q(30)*comp(19))*comp(21) 6 +q(40)*comp(11)+q(41)*comp(15)+(q(47)+q(48)+q(49))*comp(8)*comp(4) 7 +q(50)*comp(8)**2 c énergie nucléaire, on y ajoutera l'énergie graviphique epsilon(1)=SUM(epsilon(2:4)) c dérivées IF(deriv)THEN IF(h2 > 0.d0)THEN dh2t=h2*(drt(1)/r(1)-drt(2)/r(2)) dh2ro=h2*(dro(1)/r(1)-dro(2)/r(2)) DO i=1,nchim dh2x(i)=h2*(drx(1,i)/r(1)-drx(2,i)/r(2)) ENDDO dh2x(1)=dh2x(1)+h2/comp(1) ELSE dh2t=0.d0 ; dh2ro=0.d0 DO i=1,nchim dh2x(i)=0.d0 ENDDO ENDIF IF(be7 > 0.d0)THEN dent= drt(6)*mue+drt(7)*comp(1) denro=dro(6)*mue+dro(7)*comp(1) DO i=1,nchim denx(i)=drx(6,i)*mue+r(6)*dmuex(i)+drx(7,i)*comp(1) ENDDO denx(1)=denx(1)+r(7) dbe7t= be7*(drt(4)/r(4)- dent/den) dbe7ro=be7*(dro(4)/r(4)-denro/den) DO i=1,nchim dbe7x(i)=be7*(drx(4,i)/r(4)-denx(i)/den) ENDDO dbe7x(2)=dbe7x(2)+be7/comp(2); dbe7x(3)=dbe7x(3)+be7/comp(3) ELSE dbe7t=0.d0 ; dbe7ro=0.d0 DO i=1,nchim dbe7x(i)=0.d0 ENDDO ENDIF IF(li7 > 0.d0)THEN dli7t= li7*(drt(6)/r(6) +dbe7t/be7-drt(5)/r(5)) dli7ro=li7*(dro(6)/r(6)+dbe7ro/be7-dro(5)/r(5)) DO i=1,nchim dli7x(i)=li7*(drx(6,i)/r(6)+dbe7x(i)/be7 1 +dmuex(i)/mue-drx(5,i)/r(5)) ENDDO dli7x(1)=dli7x(1)-li7/comp(1) ELSE dli7t=0.d0 ; dli7ro=0.d0 DO i=1,nchim dli7x(i)=0.d0 ENDDO ENDIF c dérivées /T et=(dqt(1)*comp(1)+dqt(2)*h2+q(2)*dh2t+dqt(5)*li7+q(5)*dli7t+dqt(7)*be7 1 +q(7)*dbe7t+dqt(8)*comp(4)+dqt(9)*comp(5)+dqt(10)*comp(6)+ 2 (dqt(11)+dqt(12))*comp(7)+dqt(13)*comp(8)+dqt(14)*comp(9) 3 +dqt(34)*comp(10)+dqt(35)*comp(9)+(dqt(36)+dqt(37))*comp(18) 4 +dqt(38)*comp(14)+dqt(39)*comp(13))*comp(1) 5 +(dqt(44)+dqt(45)+dqt(46)+dqt(51))*comp(4)**2 5 +(dqt(3)*comp(2)+dqt(4)*comp(3))*comp(2)+dqt(6)*mue*be7+q(6)*mue*dbe7t 6 +(dqt(15)*comp(3)**2+dqt(16)*comp(4)+dqt(17)*comp(8)+dqt(18)*comp(6) 7 +dqt(19)*comp(10)+(dqt(20)+dqt(21))*comp(13)+dqt(25)*comp(12) 8 +dqt(31)*comp(5)+dqt(32)*comp(16)+dqt(33)*comp(17)+dqt(42)*comp(11) 9 +dqt(43)*comp(15))*comp(3) 1 +(dqt(22)*comp(13)+dqt(23)*comp(14)+dqt(24)*comp(11)+dqt(26)*comp(15) 2 +dqt(27)*comp(16)+dqt(28)*comp(17)+dqt(29)*comp(18) 3 +dqt(30)*comp(19))*comp(21)+(dqt(47)+dqt(48)+dqt(49))*comp(8)*comp(4) 4 +dqt(40)*comp(11)+dqt(41)*comp(15)+dqt(50)*comp(8)**2 c dérivées/ro ero=(dqo(1)*comp(1)+dqo(2)*h2+q(2)*dh2ro+dqo(5)*li7+q(5)*dli7ro+dqo(7)*be7 1 +q(7)*dbe7ro+dqo(8)*comp(4)+dqo(9)*comp(5)+dqo(10)*comp(6) 2 +(dqo(11)+dqo(12))*comp(7)+dqo(13)*comp(8)+dqo(14)*comp(9) 3 +dqo(34)*comp(10)+dqo(35)*comp(9)+(dqo(36)+dqo(36))*comp(18) 4 +dqo(38)*comp(14)+dqo(39)*comp(13))*comp(1) 5 +(dqo(44)+dqo(45)+dqo(46)+dqo(51))*comp(4)**2 5 +(dqo(3)*comp(2)+dqo(4)*comp(3))*comp(2) 6 +dqo(6)*mue*be7+q(6)*mue*dbe7ro+(dqo(15)*comp(3)**2+dqo(16)*comp(4) 7 +dqo(17)*comp(8)+dqo(18)*comp(6)+dqo(19)*comp(10)+(dqo(20) 8 +dqo(21))*comp(13)+dqo(25)*comp(12)+dqo(31)*comp(5)+dqo(32)*comp(16) 9 +dqo(33)*comp(17)+dqo(42)*comp(11)+dqo(43)*comp(15))*comp(3) 1 +(dqo(22)*comp(13)+dqo(23)*comp(14)+dqo(24)*comp(11)+dqo(26)*comp(15) 2 +dqo(27)*comp(16)+dqo(28)*comp(17)+dqo(29)*comp(18) 3 +dqo(30)*comp(19))*comp(21)+dqo(40)*comp(11)+dqo(41)*comp(15) 4 +(dqo(47)+dqo(48)+dqo(49))*comp(8)*comp(4)+dqo(50)*comp(8)**2 c dérivées /x(i) ex(1)=2.d0*q(1)*comp(1)+q(2)*h2+q(5)*li7+q(7)*be7 1 +q(8)*comp(4)+q(9)*comp(5)+q(10)*comp(6)+(q(11)+q(12))*comp(7) 2 +q(13)*comp(8)+q(14)*comp(9)+q(34)*comp(10)+q(35)*comp(9) 3 +(q(36)+q(37))*comp(18)+q(38)*comp(14)+q(39)*comp(13) !dérivée /H1 ex(2)=2.d0*q(3)*comp(2)+q(4)*comp(3)+q(43)*comp(15) !dérivée /He3 ex(3)=q(4)*comp(2)+3.d0*q(15)*comp(3)**2+q(16)*comp(4)+q(17)*comp(8) 1 +q(18)*comp(6)+q(19)*comp(10)+(q(20)+q(21))*comp(13)+q(25)*comp(12) 2 +q(31)*comp(5)+q(32)*comp(16)+q(33)*comp(17) 3 +q(42)*comp(11) !dérivée /He4 ex(4)=q(8)*comp(1)+q(16)*comp(3)+(q(47)+q(48)+q(49))*comp(8) 1 +2.d0*(q(44)+q(45)+q(46)+q(51))*comp(4) !dérivée /C12 ex(5)=q(9)*comp(1)+q(31)*comp(3) !dérivée /C13 ex(6)=q(10)*comp(1)+q(18)*comp(3) !dérivée /N14 ex(7)=(q(11)+q(12))*comp(1) !dérivée /N15 ex(8)=q(13)*comp(1)+q(17)*comp(3)+(q(47)+q(48)+q(49))*comp(4) 1 +2.d0*q(50)*comp(8) !dérivée /O16 ex(9)=q(14)*comp(1)+q(35)*comp(1) !dérivée /O17 ex(10)=q(19)*comp(3)+q(34)*comp(1) !dérivée /O18 ex(11)=q(24)*comp(21)+q(40)+q(42)*comp(3) !dérivée /Ne20 ex(12)=q(25)*comp(21) !dérivée /Ne21 ex(13)=q(21)*comp(3)+q(22)*comp(21)+q(39)*comp(13) !dérivée /Ne22 ex(14)=q(23)*comp(21)+q(38)*comp(1) !dérivée /Na23 ex(15)=q(26)*comp(21)+q(43)*comp(3) !dérivée /Mg24 ex(16)=q(32)*comp(3)+q(27)*comp(21)+q(41) !dérivée /Mg25 ex(17)=q(33)*comp(3)+q(28)*comp(21) !dérivée /Mg26 ex(18)=q(29)*comp(21)+(q(36)+q(37))*comp(1) !dérivée /Al27 ex(19)=q(30)*comp(21) !dérivée /Si28 ex(20)=0.d0 !dérivée /Si29 ex(21)=q(22)*comp(13)+q(23)*comp(14)+q(24)*comp(11)+q(26)*comp(15) 1 +q(27)*comp(16)+q(28)*comp(17)+q(29)*comp(18) 2 +q(30)*comp(19) !dérivée /n DO i=1,nchim !contributions des écrans ex(i)=ex(i)+(dqx(1,i)*comp(1)+dqx(2,i)*h2+q(2)*dh2x(i)+dqx(5,i)*li7 1 +q(5)*dli7x(i)+dqx(7,i)*be7+q(7)*dbe7x(i)+dqx(8,i)*comp(4) 2 +dqx(9,i)*comp(5)+dqx(10,i)*comp(6)+(dqx(11,i)+dqx(12,i))*comp(7) 3 +dqx(13,i)*comp(8)+dqx(14,i)*comp(9)+dqx(34,i)*comp(10) 4 +dqx(35,i)*comp(9)+(dqx(36,i)+dqx(37,i))*comp(18) 4 +dqx(38,i)*comp(14)+dqx(39,i)*comp(13))*comp(1) 4 +(dqx(44,i)+dqx(45,i)+dqx(46,i)+dqx(51,i))*comp(4)**2 4 +(dqx(3,i)*comp(2)+dqx(4,i)*comp(3))*comp(2)+dqx(6,i)*mue*be7 5 +q(6)*dmuex(i)*be7+q(6)*mue*dbe7x(i)+(dqx(15,i)*comp(3)**2 6 +dqx(16,i)*comp(4)+dqx(17,i)*comp(8)+dqx(18,i)*comp(6) 7 +dqx(19,i)*comp(10)+(dqx(20,i)+dqx(21,i))*comp(13)+dqx(25,i)*comp(12) 8 +dqx(31,i)*comp(5)+dqx(32,i)*comp(16)+dqx(33,i)*comp(17))*comp(3) 9 +(dqx(22,i)*comp(13)+dqx(23,i)*comp(14)+dqx(24,i)*comp(11) 1 +dqx(26,i)*comp(15)+dqx(27,i)*comp(16)+dqx(28,i)*comp(17) 2 +dqx(29,i)*comp(18)+dqx(30,i)*comp(19))*comp(21)+dqx(40,i)*comp(11) 3 +dqx(41,i)*comp(15)+dqx(42,i)*comp(11)*comp(3) 4 +dqx(43,i)*comp(15)*comp(3)+(dqx(47,i) 5 +dqx(48,i)+dqx(49,i))*comp(8)*comp(4)+dqx(50,i)*comp(8)**2 ENDDO ENDIF !deriv c-------------------------------------------------------------------------- CASE(4) !taux de production des neutrinos IF(t >= t_inf)THEN CALL rq_reac(comp,t,ro,r,drt,dro,drx,q,dqt,dqo,dqx,mue,dmuex) be7=r(4)*comp(2)*comp(3)/(r(6)*mue+r(7)*comp(1)) hhe=r(1)*comp(1)**2/amu ; be7e=r(6)*mue*be7/amu b8e=r(7)*comp(1)*be7/amu ; n13e=r(8)*comp(1)*comp(4)/amu o15e=r(10)*comp(1)*comp(6)/amu f17e=r(13)*comp(1)*comp(8)/amu ELSE hhe=0.d0 ; be7e=0.d0 ; b8e=0.d0 ; n13e=0.d0 o15e=0.d0 ; f17e=0.d0 ENDIF c-------------------------------------------------------------------- CASE DEFAULT PRINT*,'ppcno3an, fait ne peut valoir que 1, 2, 3 ou 4' PRINT*,'ERREUR fait a la valeur:',fait PRINT*,'ARRET' ; PRINT* ; STOP END SELECT RETURN END SUBROUTINE ppcno3an