Model kroplowy

Porównujemy przewidywania modelu Bethego - Weizsäckera 
                   z doświadczeniem, korzystając z IsotopeData

Formuła Bethego - Weizsäckera

Ten bardzo prosty, fenomenologiczny model jądra opiera się na obserwacji, że energia wiązania jądra jest w przybliżeniu proporcjonalna do liczby nukleonów. Nasuwa to analogię z kroplą cieczy, co z kolei wskazuje na rolę powierzchni. Trzeba też wziąć pod uwagę, że “ciecz” jądrowa jest naładowana elektrycznie, i że składa się z dwóch rodzajów fermionów. Pozwala to określić masę (energię wiązania) jądra.

Pierwszą wersję tego modelu podał w 1935 roku C. F. v. Weizsäcker. Później ulepszony i do dziś bardzo popularny, nosi nazwę modelu Bethego - Weizsäckera.

Masa nuklidu (obojętnego atomu) o liczbie masowej A i liczbie atomowej Z:

M(A,Z) = Z (Subscript[M, H])+(A-Z) (Subscript[m, n])+B(A,Z)/(c^2),

gdzie (Subscript[M, H]) to masa atomu wodoru, (Subscript[m, n]) to masa neutronu, a B (A, Z) to energia wiązania nuklidu.

(Subscript[M, H])=938.783 MeV/(c^2), (Subscript[m, n])=939.565 MeV/(c^2)

Formuła Bethego - Weizsäckera na energię wiązania:

B(A,Z)=(Subscript[b, V])A-(Subscript[b, S]) (A^(2/3))-(Subscript[b, C]) ((Z(Z-1))/(A^(1/3)))-(Subscript[b, sym])((((A-2Z))^2)/A),

Parametryzacja

Często spotykany zestaw parametrów:

bVBW=15.8;  (* współczynnik objętościowy [MeV] *)
bSBW=18.3;   (* współczynnik powierzchniowy [MeV] *)
bCBW=0.714;  (* współczynnik Coulombowski [MeV] *)
bsymBW=23.2; (* współczynnik symetrii [MeV] *)

Formuła na energię wiązania:

beBW[A_,Z_]:=bVBW*A-bSBW*A^(2/3)-bCBW*Z*(Z-1)/A^(1/3)-bsymBW*(A-2*Z)^2/A

I energia wiązania na jeden nukleon:

beBWperA[A_, Z_]:=beBW[A, Z]/A

beBWperA[4,2]

beBWperA[58,28]

Energia wiązania na nukleon

Energia wiązania na nukleon dla N = Z:

Plot[beBWperA[x,x/2],{x,2,240},
Frame->True,
FrameLabel->{"A","B/A [MeV]"},
FrameStyle->Thick,
ImageSize->700,
BaseStyle->{FontFamily->"Verdana",FontSize->22},
AxesOrigin->{0,0},
GridLines->Automatic,
PlotStyle->{{Thickness[0.005],Red}},
PlotRange->{{0,250},{0,9}}]

Obraz ogólny - obliczamy energię wiązania dla tych izotopów, które są w bazie IsotopeData :

gbBW=Table[{a-z,z, beBWperA[a,z]},{z,2,110},{a,IsotopeData[#,"MassNumber"]&/@IsotopeData[z]}]//Flatten[#,1]&;

pBW=ListPointPlot3D[gbBW,BoxRatios->{1,1,0.7},
BaseStyle->{FontFamily->"Verdana",FontSize->18},
AxesLabel->{"N","Z",""},
FaceGrids->All,
PlotRange->{{0,180},{0,110},{6,9}},
ImageSize->500, 
AxesStyle->Thick,
ColorFunction->Function[{x,y,z},Hue[(z-6)/2.5]],
ColorFunctionScaling->False
]

Porównajmy to z wartościami doświadczalnymi

gexp=Table[{a-z,z, IsotopeData[{z,a},"BindingEnergy"]},{z,2,110},{a,IsotopeData[#,"MassNumber"]&/@IsotopeData[z]}]//Flatten[#,1]&;

pexp=ListPointPlot3D[gexp,BoxRatios->{1,1,0.7},
BaseStyle->{FontFamily->"Verdana",FontSize->18},
AxesLabel->{"N","Z",""},
FaceGrids->All,
PlotRange->{{0,180},{0,110},{6,9}},
ImageSize->500, 
AxesStyle->Thick,
ColorFunction->Function[{x,y,z},Hue[(z-6)/2.5]],
ColorFunctionScaling->False
]

Show[GraphicsRow[{pBW,pexp}],ImageSize->1000]

Masy nuklidów stabilnych

Wybierzmy najpierw wszystkie nuklidy stabilne, wliczając te bardzo długożyciowe.

selectStable=Select[IsotopeData[],(IsotopeData[#,"Stable"]||IsotopeData[#,"HalfLife"]>Quantity[2*10^(16),"Seconds"])&];

Length[selectStable]

Lista {N, Z} wszystkich nuklidów stabilnych - w tej formie przyda się też później.

stableNZ={IsotopeData[#,"NeutronNumber"],
      IsotopeData[#,"AtomicNumber"]}&/@selectStable;

Aby jeszcze trochę zawęzić liczbę nuklidów, wybierzmy dla każdego pierwiastka izotop najbardziej rozpowszechniony.

Najpierw funkcja wybierająca najczęściej występujący izotop, a dokładniej liczbę masową tego izotopu :

mostAbundIzo[z_]:=Module[{iso,abu,maxabu,pos},
iso=IsotopeData[z];
abu=IsotopeData[#,"IsotopeAbundance"]&/@iso;
maxabu=Max[abu];
pos=Flatten[Position[abu,maxabu]];
IsotopeData[iso[[pos]][[1]],"MassNumber"]
]

mostAbundIzo[92]

Tworzymy listę liczb atomowych nuklidów stabilnych, bez powtórzeń:

stableZ=Union[#[[2]]&/@stableNZ];

selectMostAbund={mostAbundIzo[#],#}&/@stableZ;

Lista wartości zmierzonych {A, Bexp/A} dla tych nuklidow :

stabABperAExp={#[[1]],
      IsotopeData[{#[[2]],#[[1]]},"BindingEnergy"]}&/@selectMostAbund;

Przwidywanie formuły B-W dla tych nuklidów:

stabABperABW={#[[1]],beBWperA[#[[1]],#[[2]]]}&/@selectMostAbund;

Energia wiązania nuklidów stabilnych

Rysujemy energię wiązania na nukleon dla nuklidów stabilnych.

gstabBperAExp=ListPlot[stabABperAExp,
Frame->True,
BaseStyle->{FontFamily->"Verdana",FontSize->20},
FrameLabel->{"A","B/A [MeV]"},
PlotRange->{{0,250},{0,10}},
PlotStyle->{PointSize[0.008],Black},
FrameStyle->Thick,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
ImageSize->800
]

gstabBperABW=ListPlot[stabABperABW,
Frame->True,
BaseStyle->{FontFamily->"Verdana",FontSize->20},
FrameLabel->{"A","B/A [MeV]"},
PlotRange->{{0,250},{0,10}},
PlotStyle->{PointSize[0.01],Red},
FrameStyle->Thick,
ImageSize->800,
Joined->True]

Show[{gstabBperAExp,gstabBperABW}]

Ścieżka stabilności

Model kroplowy pozwala przewidzieć ścieżkę trwałości beta. Dla każdej liczby masowej A trzeba po prostu znaleźć izobar o najmniejszej masie.

Pełna masa nuklidu:

mBW[A_,Z_]:=Z*938.783+(A-Z)*939.565-beBW[A,Z]

Szukamy minimum tej funkcji względem Z dla ustalonego A. Obliczamy pochodną masy względem Z i badamy kiedy jest ona równa zeru:

dmdZ=D[mBW[A,Z],Z]

stZ=Solve[dmdZ==0,Z][[1]]

stableZBW[A_]:=((93.58200000000004` +((0.714`)/(A^(1/3))))/(((185.6`)/A)+((1.428`)/(A^(1/3)))))

stableZBW[208]

Możemy teraz narysować przewidywaną ścieżkę trwałości na mapie nuklidów i porównać z doświadczeniem!

tabBW=Table[{a-stableZBW[a],stableZBW[a]},{a,4,240,1}];

plotBW=ListPlot[tabBW,
Frame->True,
Joined->True,
FrameLabel->{"N","Z"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
AxesOrigin->{0,0},
PlotStyle->{{Thickness[0.004],Hue[1]}},
PlotRange->{{0,150},{0,100}},
FrameStyle->Thick,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
ImageSize->700]

Mamy już listę {N,Z} nuklidów stabilnych, zobaczmy je na mapie:

plotchart=ListPlot[stableNZ,
Frame->True,
FrameLabel->{"N","Z"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
AxesOrigin->{0,0},
PlotMarkers->{Graphics[{Black,Rectangle[{0,1}]}],0.013},
PlotRange->{{0,150},{0,100}},
FrameStyle->Thick,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
ImageSize->700]

A teraz porównajmy oba wykresy:

Show[plotBW,plotchart,ImageSize->900]

Manipulacja: analiza parametrów modelu

Zobaczmy jeszcze jak poszczególne parametry modelu B-W wpływają na energię wiązania na nukleon. Dla każdej liczby masowej wybieramy nuklid stabilny wg. modelu kroplowego i porównujemy jego zmierzoną energię wiazania z przewidywaniem.

Manipulate[
quality=0;
BWdata=Table[{a,be[a,Round[stabZ[a]]]},{a,5,239,2}];
(* Do[quality=quality+(expdata[[i,2]]-BWdata[[i,2]])^2,{i,1,118,1}];*)
quality=Total[(expdata[[All,2]]-BWdata[[All,2]])^2];
qq=Round[quality/2.49006,0.001];
gexp=ListPlot[expdata,PlotStyle->{PointSize[0.008],Black}];
gVol=Plot[bV,{x,1,240},
PlotStyle->{{Thickness[0.005],Hue[0.55]}}];
gVolS=Plot[bV-bS/x^(1/3),{x,1,240},
PlotStyle->{{Thickness[0.005],Hue[0.7]}}];
gVolSC=Plot[bV-bS/x^(1/3)-bC*stabZ[x]*(stabZ[x]-1)/x^(4/3),{x,1,240},
PlotStyle->{{Thickness[0.005],Hue[0.8]}}];
gfull=Plot[be[x,stabZ[x]],{x,1,240},
PlotStyle->{{Thickness[0.006],Hue[1]}}];
Show[gVol,gVolS,gVolSC,gfull,gexp,
Frame->True,
AxesOrigin->{0,0},
PlotRange->{{1,240},{4,17}},
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
FrameLabel->{"A","B/A [MeV]"},
BaseStyle->{FontFamily->"Verdana",FontSize->22},
ImageSize->750,
Epilog->{Text[Style[StringForm["Q = `1`",qq],Red,24],{115,5}],Text[Style["V",Hue[0.55],22],{220,16.3}],
Text[Style["V+S",Hue[0.7],22],{220,13.5}],
Text[Style["V+S+C",Hue[0.8],22],{220,10}],
Text[Style["V+S+C+sym",Hue[1],22],{220,6}]
}],
{{bV,15.8,Style["(Subscript[b, V]) ",Hue[0.55],24]},10,20,0.1,Appearance->"Labeled",LabelStyle->Directive[Hue[0.55],24]},
{{bS,18.3,Style["(Subscript[b, S]) ",Hue[0.7],24]},10,25,0.1,Appearance->"Labeled", LabelStyle->Directive[Hue[0.7],24]},
{{bC,0.714,Style["(Subscript[b, C]) ",Hue[0.8],24]},0.3,1.5,0.001,Appearance->"Labeled",
LabelStyle->Directive[Hue[0.8],24]},
{{bsym,23.2,Style["(Subscript[b, sym])",Hue[1.0],24]},10,40,0.1,Appearance->"Labeled",
LabelStyle->Directive[Hue[1],24]},
TrackedSymbols->{bV,bS,bC,bsym},
FrameLabel->{None, None, fulltit,None},
SaveDefinitions->True,
Initialization:>{
be[A_,Z_]:=(bV*A-bS*A^(2/3)-bC*Z*(Z-1)/A^(1/3)-bsym*(A-2*Z)^2/A)/A;
stabZ[A_]:=((93.58200000000004` +((0.714`)/(A^(1/3))))/(((185.6`)/A)+((1.428`)/(A^(1/3)))));
expdata=Table[{a, QuantityMagnitude[IsotopeData[{Round[stabZ[a]],a},"BindingEnergy"]]},{a,5,239,2}];titleA=Text[Style["Energia wiązania na nukleon w modelu kroplowym",Blue,FontFamily->"Verdana",24]];
titleB=Text[Style["B(A,Z)/A = (Subscript[b, V]) - (Subscript[b, S])(1/(A^(1/3))) - (Subscript[b, C])((Z (Z - 1))/(A^(4/3))) - (Subscript[b, sym])((((A - 2  Z))^2)/(A^2))",FontFamily->"Georgia",Purple,24]];
fulltit=GraphicsColumn[{titleA,titleB},AspectRatio->0.15,ImageSize->700,Frame->True]}
]

Masy wszystkich nuklidów

Spoglądamy teraz na całą mapę nuklidów. Żeby sprawdzić model kroplowy obliczamy różnicę między zmierzoną energią wiązania na nukleon a przewidywaniem modelu B-W dla wszystkich nuklidów:

diff=Table[{a-z,z, QuantityMagnitude[IsotopeData[{z,a},"BindingEnergy"]]-beBWperA[a,z]},{z,1,110},{a,IsotopeData[#,"MassNumber"]&/@IsotopeData[z]}];

diff[[1]]

gdiff=ListDensityPlot[Flatten[diff,1], 
PlotRange->{{10,140},{10,100},{-0.15,0.15}},
ColorFunction->ColorData["SunsetColors"],
ColorFunctionScaling->True,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
Frame->True,
FrameLabel->{{"Z",None},{"N","B/(Subscript[A, exp])-B/(Subscript[A, BW])"}},
BaseStyle->{FontFamily->"Verdana",FontSize->22},
AspectRatio->0.75,
PlotLegends->Placed[Automatic,Right],
ImageSize->700]

Efekt sił pairing

Rozważmy obszar wokół izotopów cyny:

gdiff2=ListDensityPlot[Flatten[diff,1], 
PlotRange->{{49,85},{40,60},{-0.1,0.12}},
ColorFunction->ColorData["SunsetColors"],
ColorFunctionScaling->True,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
Frame->True,
FrameLabel->{{"Z",None},{"N","B/(Subscript[A, exp])-B/(Subscript[A, BW])"}},
BaseStyle->{FontFamily->"Verdana",FontSize->22},
AspectRatio->0.7,
PlotLegends->Placed[Automatic,Right],
ImageSize->600]

Zobaczmy dokładniej różnice między eksperymentem a modelem dla izotopów cyny (Z=50) i antymonu (Z=51):

diffZ50={#[[1]],#[[3]]}&/@diff[[50]];
diffZ51={#[[1]],#[[3]]}&/@diff[[51]];

g50=ListPlot[diffZ50,
Frame->True,
Joined->True,
MeshStyle->Red,
Mesh->All,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
FrameLabel->{"N","(B/A)_exp-(B/A)_BW [MeV]"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotRange->{Automatic,{-0.05,0.2}},
PlotStyle->GrayLevel[0.5],
ImageSize->700];

g51=ListPlot[diffZ51,
Frame->True,
Joined->True,
MeshStyle->Blue,
Mesh->All,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
FrameLabel->{"N","(B/A)_exp-(B/A)_BW [MeV]"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotRange->{Automatic,{-0.05,0.2}},
PlotStyle->GrayLevel[0.5],
ImageSize->700];

Show[g50,g51]

Wyznaczamy (Subscript[Δ, n])

Dla nieparzystego neutronu:

deltaN[n_,z_]:=Module[{x0,xm1,xp1,a,d},
If[EvenQ[n]&&EvenQ[z],
a=n+z;
x0=QuantityMagnitude[IsotopeData[{z,a},"BindingEnergy"]]*a;
xm1=QuantityMagnitude[IsotopeData[{z,a-1},"BindingEnergy"]]*(a-1);
xp1=QuantityMagnitude[IsotopeData[{z,a+1},"BindingEnergy"]]*(a+1);
If[NumericQ[x0]&&NumericQ[xm1]&&NumericQ[xp1],
d=(2*x0-xm1-xp1)/2,
d=0],
d=0];
d
]

deltaN[20,20]

Dla izotonów zadanego N:

fdeltaN[nn_]:=Module[{lizo,m,lfin},
lizo=Select[IsotopeData[],(IsotopeData[#,"NeutronNumber"]==nn)&& (EvenQ[IsotopeData[#,"AtomicNumber"]])&];
      lfin={nn,deltaN[nn,IsotopeData[#,"AtomicNumber"]]}&/@lizo;
      Select[lfin,#[[2]]>0&]
]

Zbieramy dane dla wszystkich parzystych N i rysujemy:

tN=Table[fdeltaN[n],{n,10,150,2}];

pldN=ListPlot[tN,
Joined->True,
Mesh->All,
MeshStyle->Directive[Darker[Blue],PointSize[0.007]],
FrameLabel->{"N","(Subscript[Δ, n]) [MeV]"},
Frame->True,
PlotRange->{{5,152},{0,5}},                      BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotStyle->Automatic,
ImageSize->700,
Filling->False
              ]

Jak wiadomo, jako wartość Δ przyjmuje się Δ=(12/(Sqrt[A]))MeV. Zobaczmy jak to się ma do danych doświadczalnych. Wartość Z bierzemy z funkcji stableZBW[A] :

fdN=Table[{a-Round[stableZBW[a]],12/Sqrt[a]},{a,10,250,2}];

plfdN=ListPlot[fdN,
Joined->True,
FrameLabel->{"N","(Subscript[Δ, n]) [MeV]"},
Frame->True,
PlotRange->{{5,150},{0,5}},    BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotStyle->{{Thickness[0.005],Hue[1]}},
ImageSize->700,
Filling->False];

Show[{pldN,plfdN}]

Wyznaczamy (Subscript[Δ, p])

Powtarzamy te same kroki dla protonów.
Dla nieparzystego protonu:

deltaZ[n_,z_]:=Module[{x0,xm1,xp1,a,d},
If[EvenQ[n]&&EvenQ[z],
a=n+z;
x0=QuantityMagnitude[IsotopeData[{z,a},"BindingEnergy"]]*a;
xm1=QuantityMagnitude[IsotopeData[{z-1,a-1},"BindingEnergy"]]*(a-1);
xp1=QuantityMagnitude[IsotopeData[{z+1,a+1},"BindingEnergy"]]*(a+1);
If[NumericQ[x0]&&NumericQ[xm1]&&NumericQ[xp1],
d=(2*x0-xm1-xp1)/2,
d=0],
d=0];
d
]

deltaZ[60,40]

Dla izotopów danego Z :

fdeltaZ[zz_]:=Module[{lizo,m,lfin},
lizo=Select[IsotopeData[],(IsotopeData[#,"AtomicNumber"]==zz)&& (EvenQ[IsotopeData[#,"NeutronNumber"]])&];
      lfin={zz,deltaZ[IsotopeData[#,"NeutronNumber"],zz]}&/@lizo;
      Select[lfin,#[[2]]>0&]
]

|

fdeltaZ[30]

Dla wszystkich Z:

tZ=Table[fdeltaZ[z],{z,10,92,2}];

pldZ=ListPlot[tZ,
Joined->True,
Mesh->All,
MeshStyle->Directive[Darker[Blue],PointSize[0.007]],
FrameLabel->{"Z","(Subscript[Δ, p]) [MeV]"},
Frame->True,
PlotRange->{{5,94},{0,5}},                      BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotStyle->Automatic,
ImageSize->700,
Filling->False
              ]

fdZ=Table[{stableZBW[a],16/Sqrt[a]},{a,10,240,2}];

plfdZ=ListPlot[fdZ,
Joined->True,
FrameLabel->{"Z","(Subscript[Δ, p]) [MeV]"},
Frame->True,
PlotRange->{{5,94},{0,5}},    BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotStyle->{{Thickness[0.005],Hue[1]}},
ImageSize->700,
Filling->False];

Show[{pldZ,plfdZ}]

Pełny model z członem pairing

Dodajemy człon opisujący siły pairing do modelu B-W.

bpairBW=12.0; (* parametr członu pairing [MeV] *)

beBWpair[A_,Z_]:=beBW[A,Z]+If[EvenQ[A],If[EvenQ[Z],1 ,-1 ],0]*bpairBW/Sqrt[A]

Sprawdźmy jak teraz wygląda różnica między doświadczeniem a modelem B-W dla energii wiązania:

diffpair=Table[{a-z,z, QuantityMagnitude[IsotopeData[{z,a},"BindingEnergy"]]-beBWpair[a,z]/a},{z,1,110},{a,IsotopeData[#,"MassNumber"]&/@IsotopeData[z]}];

gdiffpair=ListDensityPlot[Flatten[diffpair,1], 
PlotRange->{{10,140},{10,100},{-0.15,0.15}},
ColorFunction->ColorData["SunsetColors"],
ColorFunctionScaling->True,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
Frame->True,
FrameLabel->{{"Z",None},{"N","B/(Subscript[A, exp])-B/(Subscript[A, BW])"}},
BaseStyle->{FontFamily->"Verdana",FontSize->22},
AspectRatio->0.75,
PlotLegends->Placed[Automatic,Right],
ImageSize->700]

Obszar wokół cyny

Wybieramy jeszcze raz ten sam obszar wokół cyny:

gdiff2pair=ListDensityPlot[Flatten[diffpair,1], 
PlotRange->{{49,85},{40,60},{-0.1,0.12}},
ColorFunction->ColorData["SunsetColors"],
ColorFunctionScaling->True,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
Frame->True,
FrameLabel->{{"Z",None},{"N","B/(Subscript[A, exp])-(B/(Subscript[A, BW])+pairing)"}},
BaseStyle->{FontFamily->"Verdana",FontSize->22},
AspectRatio->0.7,
PlotLegends->Placed[Automatic,Right],
ImageSize->600];

Show[GraphicsRow[{gdiff2,gdiff2pair}],ImageSize->1100]

Dla samych izotopów cyny i antymonu :

diffZ50pair={#[[1]],#[[3]]}&/@diffpair[[50]];
diffZ51pair={#[[1]],#[[3]]}&/@diffpair[[51]];

g50pair=ListPlot[diffZ50pair,
Frame->True,
Joined->True,
MeshStyle->Red,
Mesh->All,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
FrameLabel->{"N","(B/A)_exp-(B/A)_(BW + pair) [MeV]"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotRange->{Automatic,{-0.05,0.2}},
PlotStyle->GrayLevel[0.5],
ImageSize->700];

g51pair=ListPlot[diffZ51pair,
Frame->True,
Joined->True,
MeshStyle->Blue,
Mesh->All,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
FrameLabel->{"N","(B/A)_exp-(B/A)_(BW + pair) [MeV]"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotRange->{Automatic,{-0.05,0.2}},
PlotStyle->GrayLevel[0.5],
ImageSize->700];

Show[GraphicsColumn[{Show[g50,g51],Show[g50pair,g51pair]}],ImageSize->600]

Manipulacja : porównanie dla dowolnego Z

Manipulate[
dataZ=diff[[z]];
ListPlot[dataZ,
Frame->True,
Joined->True,
MeshStyle->Red,
Mesh->All,
GridLines->Automatic,
GridLinesStyle->Directive[Dotted],
FrameLabel->{"N","(B/A)_exp-(B/A)_(BW + pair) [MeV]"},
BaseStyle->{FontFamily->"Verdana",FontSize->20},
PlotRange->Automatic,
PlotStyle->GrayLevel[0.5],
ImageSize->700
],
{{z,50,Style["Z",Purple,22]},10,110,1,Appearance->"Labeled",LabelStyle->Directive[Purple,22]},
TrackedSymbols->{z},
FrameLabel->{None, None, fulltit,None},
SaveDefinitions->True,
Initialization:>{
bepair[a_,z_]:=(15.8*a-18.3*a^(2/3)-0.714*z*(z-1)/a^(1/3)-23.2*(a-2*z)^2/a +If[EvenQ[a],If[EvenQ[z],1 ,-1 ],0]*12.0/Sqrt[a])/a;
diff=Table[{a-z, QuantityMagnitude[IsotopeData[{z,a},"BindingEnergy"]]-bepair[a,z]},{z,1,110},{a,IsotopeData[#,"MassNumber"]&/@IsotopeData[z]}];
titleA=Text[Style["Energia wiązania na nukleon",Blue,FontFamily->"Verdana",24]];
titleB=Text[Style["Porównanie modelu kroplowego z eksperymentem",FontFamily->"Verdana",Black,22]];
fulltit=GraphicsColumn[{titleA,titleB},AspectRatio->0.10,ImageSize->700,Frame->True]},
SynchronousInitialization->False
]