(* Informację o tempach reakcji termojądrowym można pobrać np: z http://ie.lbl.gov/astro/friedel.html *)

rate33txt=Import["https://ie.lbl.gov/astro2/rate6.txt","Lines"][[15]]

ratepeptxt=Import["https://ie.lbl.gov/astro2/rate4.txt","Lines"][[13]]

ratepptxt=Import["https://ie.lbl.gov/astro2/rate4.txt","Lines"][[14]]

ratepdtxt=Import["https://ie.lbl.gov/astro2/rate4.txt","Lines"][[16]]

(* Funkcja importująca text wyprodukowany przez prawdopodobnie FOTRAN77 *)
importRate[rate_?StringQ]:=
{
ImportString[StringTake[rate,{-92,-80}],"Table"][[1,1]],
ImportString[StringTake[rate,{-79,-67}],"Table"][[1,1]],
ImportString[StringTake[rate,{-66,-54}],"Table"][[1,1]],
ImportString[StringTake[rate,{-53,-40}],"Table"][[1,1]],
ImportString[StringTake[rate,{-39,-27}],"Table"][[1,1]],
ImportString[StringTake[rate,{-26,-14}],"Table"][[1,1]],
ImportString[StringTake[rate,{-13,-1 } ],"Table"][[1,1]]
}

(* Standardowy sposób interpolacji w astrofizyce *)
rate[a_?ListQ]:=Exp[a[[1]]+a[[2]]/T9+a[[3]]/T9^(1/3)+a[[4]] T9^(1/3)+a[[5]] T9+a[[6]] T9^(5/3)+a[[7]] Log[T9]]

rate[{a1,a2,a3,a4,a5,a6,a7}] (* Struktura wzoru, T9 = T/((10^9) K) *)

(* Układ równań różniczkowych opisujących spalanie wodoru w cyklu ppI *)

eq1=np'[t]==-2 λpp np[t]^2 - λpd np[t] nd[t]+2 λ33 n3[t]^2

eq2=nd'[t]==λpp np[t]^2 - λpd np[t] nd[t]

eq3=n3'[t]==λpd np[t] nd[t]-2 λ33 n3[t]^2

eq4=nα'[t]==λ33 n3[t]^2

(* Spr: zachowanie liczby barionowej *)
eq1[[2]]+2 eq2[[2]]+3eq3[[2]]+4 eq4[[2]]//Simplify

n={np[t],nd[t],n3[t],nα[t]}

f={eq1[[2]],2 eq2[[2]],3eq3[[2]],4 eq4[[2]]}

jac=Table[D[f[[i]],n[[j]]],{i,1,4},{j,1,4}]

NA=6.022*10^23 (* Liczba Avogadro *)

(* Wzory z podręcznika Bisnovatyj-Kogan, Stellar Physics, Tom I, równanie (4.2.5) na str. 139 (4.2 Hydrogen, Deuterium and Helium Burning ) *)

(*pp=4.01*10^-15 T9^(-2/3) Exp[-3.380/T9^(1/3)]*(1+0.123 T9^(1/3)+1.09*T9^(2/3)+0.938 T9)*)

(* Wzory z bazy danych ie.lbl.gov *)
λpp=rate[importRate[ratepptxt]]/NA

(*pd=2.24*10^3T9^(-2/3) Exp[-3.72/T9^(1/3)]*(1+0.112 T9^(1/3)+3.38*T9^(2/3)+2.65 T9)*)

λpd=rate[importRate[ratepdtxt]]/NA

(*λ33=6.04*10^10T9^(-2/3) Exp[-12.276/T9^(1/3)]*(1+0.034 T9^(1/3)-0.522*T9^(2/3)-0.124 T9+0.353 T9^(4/3)+0.213 T9^(5/3))*)

λ33=rate[importRate[rate33txt]]/NA

LogLogPlot[{λpp,λpd,λ33},{T9,0.001,0.1},PlotStyle->{Red,Blue,Green}]

T9=0.015 (* Temperatura w centrum Słońca, 15 milionów K *)

ρ=150.0 (* Gęstość w centrum Słońca *)

endTime=10^11 * 365*86400. (* Czas spalania: 100 miliardów lat *)

{λpp,λpd,λ33}

(* Rozwiązanie układu równań różniczkowych. Proszę zwrócić uwagę na użycie metody typu implicit "BDF" oraz podniesienie precyzji obliczeń do 32 cyfr znaczących *)
rozw=NDSolve[{eq1,eq2,eq3,eq4,np[0]==12/13*NA*ρ,nd[0]==0.0 NA ρ,n3[0]==0,nα[0]==1/13NA*ρ},{np,nd,n3,nα},{t,0,endTime},MaxSteps->100000,Method->"BDF",WorkingPrecision->32][[1]]

npF=np/.rozw;
ndF=nd/.rozw;
n3F=n3/.rozw;
nαF=nα/.rozw;

mldLat=10^9 * 365*86400.0 (* 1 miliard lat w sekundach *)

w1=Plot[{Tooltip[npF[t mldLat],"H"],Tooltip[ndF[t mldLat],"D"],Tooltip[n3F[t mldLat],"He3"],Tooltip[nαF[t mldLat],"He4"]},{t,0,10},PlotRange->All,PlotLegends->{"Wodór", "Deuter" ,"Hel 3","Hel 4" },AxesLabel->{"Czas [mld lat]","Ilość jąder atomów [1/(cm^3)]"},BaseStyle->FontSize->14,ImageSize->600]

Export["SpalanieH_ppI_1.pdf",w1]

w2=LogLogPlot[{Tooltip[npF[t],"H"],Tooltip[ndF[t],"D"],Tooltip[n3F[t],"He3"],Tooltip[nαF[t],"He4"]},{t,10^-22,endTime},PlotRange->All,PlotLegends->{"Wodór", "Deuter" ,"Hel 3","Hel 4" },AxesLabel->{"Czas","Ilość jąder atomów [1/(cm^3)]"},BaseStyle->FontSize->14,ImageSize->600,Ticks->{{{1,"1 s"},{86400,"1\ndzień"},{1000*365.25*86400,"1000\nlat"},{mldLat,"1\nmld lat"}},Automatic}]

Export["SpalanieH_ppI_2.pdf",w2]

Xp = npF[t]/(npF[t]+ndF[t]+n3F[t]+nαF[t]);

Xd = ndF[t]/(npF[t]+ndF[t]+n3F[t]+nαF[t]);

X3 = n3F[t]/(npF[t]+ndF[t]+n3F[t]+nαF[t]);

Xα = nαF[t]/(npF[t]+ndF[t]+n3F[t]+nαF[t]);

LogLinearPlot[{Xp,Xd,X3,Xα},{t,1,endTime},PlotRange->{0,1}]

LogLogPlot[npF[t]+2ndF[t]+3n3F[t]+4nαF[t],{t,1,endTime}] (* Test zachowania liczby barionowej przez ALGORYTM NDSolve *)

Plot[{Exp[-e],Exp[-1/Sqrt[e]],Exp[-e]*Exp[-1/Sqrt[e]]},{e,0,6},PlotRange->All,PlotLegends->{"(e^(-E/kT))", "(e^(-(Subscript[Z, 1]) (Subscript[Z, 2])/√E))","Pik Gamowa"}]

Export["PikGamowa.pdf",%]

jac/.rozw/.t->endTime/2//MatrixPlot