dataAKB=Import["C:\\Users\\Natalia Morais\\MEOCloud\\problema10e11\\pp_si_vkbofk_new.dat","Table"];
maxplot=Position[dataAKB,10.`][[1]][[1]];
dataplot=Take[dataAKB,maxplot];
ListPlot[{dataplot},PlotRange->All,AxesLabel->{"k","A00"}]

interAKBk=Interpolation[dataAKB];
Plot[interAKBk[k],{k,0.01,10},PlotRange->All,AxesLabel->{"k","A00"}]

model[k_,B_,Rb_]:=(1/4) B (E^(-(1/4) (k^2) (Rb^2))) (Sqrt[π]) (Rb^3)

fit=FindFit[dataAKB,model[k,B,Rb],{B,Rb},k]

Show[Plot[model[k,B,Rb]/.fit,{k,0,10},PlotRange->All],ListPlot[dataplot,PlotStyle->Purple],AxesLabel->{"k","A00"}]

min[B_,Rb_]:=NIntegrate[(((interAKBk[q]-model[q,B,Rb]))^2)*(q^2),{q,0.01,15}]

Dynamic[B]
Dynamic[Rb]
minRb=0.8;
maxRb=1.3;
dRb=0.01;
minB=4;
maxB=13;
dB=0.5;
contour= Table[{Rb,B,min[B,Rb]},{Rb,minRb,maxRb,dRb},{B,minB,maxB,dB}];

countour2=Flatten[contour,1];

Show[ListContourPlot[countour2,PlotRange->{0,5},ColorFunction->"Rainbow"],ListPlot[{{Rb/.fit,B/.fit}},PlotStyle->Green],FrameLabel->{"Rb","B"}]

vldadata=Import["C:\\Users\\Natalia Morais\\MEOCloud\\problema10e11\\pp_si_vofr_new.dat","Table"];
vldas=Table[{vldadata[[i]][[1]],vldadata[[i]][[2]]+vldadata[[i]][[4]]},{i,Dimensions[vldadata][[1]]}];
intervldas=Interpolation[vldas];
vldap=Table[{vldadata[[i]][[1]],vldadata[[i]][[3]]+vldadata[[i]][[4]]},{i,Dimensions[vldadata][[1]]}];
intervldap=Interpolation[vldap];
dataAKB=Import["C:\\Users\\Natalia Morais\\MEOCloud\\problema10e11\\pp_si_vkbofr_new.dat","Table"];
interAKB=Interpolation[dataAKB];
Plot[{intervldas[x],intervldap[x],interAKB[x]},{x,0,5},AxesLabel->{"r",},PlotLegends->{"Vs","Vp","a00"}]

sol=ParametricNDSolve[{-(1/2)y''[x]+intervldas[x] y[x]==En y[x],y[0.01]==0.01,y'[0.01]==1},y,{x,0.01,20},{En}];
root=FindRoot[y[En][19.9]/.sol,{En,-0.398}]
en=En/.root;
fsnn=y[en]/.sol;
norm=Sqrt[NIntegrate[(fsnn[r])^2,{r,0.01,10}]];
fs[x_] = fsnn[x] / norm;

Plot[fs[x],{x,0.01,5},AxesLabel->{"r","Psi00"}]

norm0 =Sqrt[NIntegrate[fs[x]^2,{x,0,10}]];
f1z[x_]=fs[x] / norm0;
en1x=-1;
en1=0;

Dynamic[en1]
For[i=0,Abs[en1-en1x]>(10^(-6)),++i,
en1x=en1;
prod =NIntegrate[interAKB[r]f1z[r]r,{r,0.01,10}];
sol1=ParametricNDSolve[{-(1/2)y1''[x]+intervldap[x] y1[x]+x interAKB[x]prod==En  y1[x],y1[0.01]==f1z[0.01],y1'[0.01]==f1z'[0.01]},y1,{x,0.01,10},{En}];
root1=FindRoot[y1[En][9.9]/.sol1,{En,-0.398}];
en1=En/.root1;
f1nn=y1[en1]/.sol1;
norm1=Sqrt[NIntegrate[(f1nn[r])^2,{r,0.01,10}]];
f1[x_] = f1nn[x] / norm1;
f1x[x_]=f1z[x];
f1znn[x_]=0.7f1[x]+0.3f1x[x];
norm1=Sqrt[NIntegrate[(f1znn[r])^2,{r,0.01,10}]];
f1z[x_] = f1znn[x] / norm1;
]

Clear[Ra,ktf,qz]
siPseudo[k_,Ra_,ktf_,qz_]=16π(((k^2)/(qz^2))-1)((E^(-(((k^2)(Ra^2))/4)))/((k^2)+(ktf^2)));
Vlocal[r_,Ra_,ktf_,qz_]:=NIntegrate[siPseudo[k,Ra,ktf,qz]k^2 SphericalBesselJ[0,k r],{k,0,10}]/(2 Pi^2)

Ra=0.93;
ktf=0.527449;
qz=2.17;
tVlocal=Table[{r,Vlocal[r,Ra,ktf,qz]},{r,0,20,0.1}];
intertVlocal=Interpolation[tVlocal];

Plot[{intertVlocal[x],intervldap[x]},{x,0,5},AxesLabel->{"r","Vlocal"},PlotLegends->{"our model","from reference"}]

Clear[f2,f2z,f2x,f2nn,f2znn]
norm0 =Sqrt[NIntegrate[fs[x]^2,{x,0,10}]];
f2z[x_]=fs[x] / norm0;
en2x=-1;
en2=0;

Dynamic[en2]
For[i=0,Abs[en2-en2x]>(10^(-6)),++i,
en2x=en2;
prod =NIntegrate[interAKB[r]f2z[r]r,{r,0.01,10}];
sol2=ParametricNDSolve[{-(1/2)y2''[x]+intertVlocal[x] y2[x]+x interAKB[x]prod==En  y2[x],y2[0.01]==f2z[0.01],y2'[0.01]==f2z'[0.01]},y2,{x,0.01,10},{En}];
root2=FindRoot[y2[En][9.9]/.sol2,{En,-0.398}];
en2=En/.root2;
f2nn=y2[en2]/.sol2;
norm2=Sqrt[NIntegrate[(f2nn[r])^2,{r,0.01,10}]];
f2[x_] = f2nn[x] / norm2;
f2x[x_]=f2z[x];
f2znn[x_]=0.7f2[x]+0.3f2x[x];
norm2=Sqrt[NIntegrate[(f2znn[r])^2,{r,0.01,10}]];
f2z[x_] = f2znn[x] / norm2;
]

Clear[ktf]
minktf=0.535;
maxktf=0.545;
dktf=0.001;
Ektf=Table[,{k,minktf,maxktf,dktf}];
maxj=Dimensions[Ektf][[1]];

Dynamic[ktf]
Dynamic[en2]
Do[
ktf=minktf+(j-1)dktf;
tVlocal=Table[{r,Vlocal[r,Ra,ktf,qz]},{r,0,20,0.1}];
intertVlocal=Interpolation[tVlocal];
Clear[f2,f2z,f2x,f2nn,f2znn];
norm0 =Sqrt[NIntegrate[fs[x]^2,{x,0,10}]];
f2z[x_]=fs[x] / norm0;
en2x=-1;
en2=0;
For[i=0,Abs[en2-en2x]>(10^(-6)),++i,
en2x=en2;
prod =NIntegrate[interAKB[r]f2z[r]r,{r,0.01,10}];
sol2=ParametricNDSolve[{-(1/2)y2''[x]+intertVlocal[x] y2[x]+x interAKB[x]prod==En  y2[x],y2[0.01]==f2z[0.01],y2'[0.01]==f2z'[0.01]},y2,{x,0.01,10},{En}];
root2=FindRoot[y2[En][9.9]/.sol2,{En,-0.398}];
en2=En/.root2;
f2nn=y2[en2]/.sol2;
norm2=Sqrt[NIntegrate[(f2nn[r])^2,{r,0.01,10}]];
f2[x_] = f2nn[x] / norm2;
f2x[x_]=f2z[x];
f2znn[x_]=0.7f2[x]+0.3f2x[x];
norm2=Sqrt[NIntegrate[(f2znn[r])^2,{r,0.01,10}]];
f2z[x_] = f2znn[x] / norm2;
];

Ektf[[j]]={ktf,en1-en2}
,{j,maxj}]

ListLinePlot[Ektf]

interEktf=Interpolation[Ektf];
solktf=FindRoot[interEktf[x],{x,0.539}]
ktf=x/.solktf;

tVlocal=Table[{r,Vlocal[r,Ra,ktf,qz]},{r,0,20,0.1}];
intertVlocal=Interpolation[tVlocal];

Clear[f2,f2z,f2x,f2nn,f2znn]
norm0 =Sqrt[NIntegrate[fs[x]^2,{x,0,10}]];
f2z[x_]=fs[x] / norm0;
en2x=-1;
en2=0;

Dynamic[en2]
For[i=0,Abs[en2-en2x]>(10^(-6)),++i,
en2x=en2;
prod =NIntegrate[interAKB[r]f2z[r]r,{r,0.01,10}];
sol2=ParametricNDSolve[{-(1/2)y2''[x]+intertVlocal[x] y2[x]+x interAKB[x]prod==En  y2[x],y2[0.01]==f2z[0.01],y2'[0.01]==f2z'[0.01]},y2,{x,0.01,10},{En}];
root2=FindRoot[y2[En][9.9]/.sol2,{En,-0.398}];
en2=En/.root2;
f2nn=y2[en2]/.sol2;
norm2=Sqrt[NIntegrate[(f2nn[r])^2,{r,0.01,10}]];
f2[x_] = f2nn[x] / norm2;
f2x[x_]=f2z[x];
f2znn[x_]=0.7f2[x]+0.3f2x[x];
norm2=Sqrt[NIntegrate[(f2znn[r])^2,{r,0.01,10}]];
f2z[x_] = f2znn[x] / norm2;
]

qz
Ra
ktf

Clear[B,Rb]
model[k_,B_,Rb_]:=(1/4) B (E^(-(1/4) (k^2) (Rb^2))) (Sqrt[π]) (Rb^3)
Vnl[r_,B_,Rb_]:=NIntegrate[model[k,B,Rb]k^2 SphericalBesselJ[0,k r],{k,0,7}]/(Sqrt[Pi]Pi)

Dynamic[r]
Dynamic[B]
Dynamic[Rb]
tVnl=Table[{{r,B,Rb},Vnl[r,B,Rb]},{r,0,10,0.2},{B,minB,maxB,dB},{Rb,minRb,maxRb,dRb}];

datatVnl=Flatten[tVnl,2];
interVnl=Interpolation[datatVnl];

prod0=NIntegrate[interAKB[r]fs[r]r,{r,0.01,10}]

prod1=Table[{Rb,B,NIntegrate[interVnl[r,B,Rb]fs[r]r,{r,0.01,10}]},{Rb,minRb,maxRb,dRb},{B,minB,maxB,dB}];
maxi=Dimensions[prod1][[1]]
maxj=Dimensions[prod1][[2]]

difp=Table[{prod1[[i]][[j]][[1]],prod1[[i]][[j]][[2]],Abs[prod1[[i]][[j]][[3]]-prod0]},{i,maxi},{j,maxj}];

difpcontour=Flatten[difp,1];
ListContourPlot[difpcontour,PlotRange->{0,0.2},FrameLabel->{"Rb","B"}]

dataprod=Flatten[prod1,1];
interprod=Interpolation[dataprod]

valBRb=Table[,{Rb,minRb,maxRb,dRb}];
maxi=Dimensions[valBRb][[1]];

minBRb=Table[{Rb=minRb+(i-1)*dRb,difp[[i]][[Position[difp[[i]],Min[difp[[i]]]][[1]][[1]]]][[2]]},{i,maxi}];
interminBRb=Interpolation[minBRb]

norm0 =Sqrt[NIntegrate[fs[x]^2,{x,0,10}]];
f3z[x_]=fs[x] / norm0;

Dynamic[Rb]
Do[

Rb=minRb+(i-1)*dRb;

sol3[B_]:=NDSolve[{-(1/2)y3''[x]+intertVlocal[x] y3[x]+x interVnl[x,B,Rb]interprod[Rb,B]==-0.398 y3[x],y3[0.01]==f3z[0.01],y3'[0.01]==f3z'[0.01]},y3,{x,0.01,10}];

root3=FindRoot[y3[10] /. sol3[B],{B,interminBRb[Rb]}];

valBRb[[i]]={Rb,B/.root3};

,{i,maxi}]
Clear[Rb]

Show[ListContourPlot[difpcontour,PlotRange->{0,0.2}],ListPlot[{{Rb/.fit,B/.fit}},PlotStyle->Green],ListPlot[valBRb,PlotStyle->{Red,PointSize[0.01]}],FrameLabel->{"Rb","B"}]

valBRb

ListPlot[valBRb,PlotStyle->{Red,PointSize[0.01]},AxesLabel->{"Rb","B"}]

Clear[Rb]

fB[Rb_]=Fit[valBRb,{1,Rb,Rb^2},Rb]

Show[ListContourPlot[countour2,PlotRange->{0,5},ColorFunction->"Rainbow"],ListPlot[{{Rb/.fit,B/.fit}},PlotStyle->Green],ListPlot[valBRb,PlotStyle->{Red}],Plot[fB[Rb],{Rb,minRb,maxRb},PlotStyle->Yellow],FrameLabel->{"Rb","B"}]

valBRbM=Table[,{Rb,minRb,maxRb,dRb}];
maxi=Dimensions[valBRb][[1]];

Dynamic[Rb]
Do[

Rb=minRb+(i-1)*dRb;

sol3[B_]:=NDSolve[{-(1/2)y3''[x]+intertVlocal[x] y3[x]+x interVnl[x,B,Rb]interprod[Rb,B]==(-0.398+0.5/27.211) y3[x],y3[0.01]==f3z[0.01],y3'[0.01]==f3z'[0.01]},y3,{x,0.01,10}];

root3=FindRoot[y3[10] /. sol3[B],{B,interminBRb[Rb]}];

valBRbM[[i]]={Rb,B/.root3};

,{i,maxi}]

valBRbm=Table[,{Rb,minRb,maxRb,dRb}];
maxi=Dimensions[valBRb][[1]];

Dynamic[Rb]
Do[

Rb=minRb+(i-1)*dRb;

sol3[B_]:=NDSolve[{-(1/2)y3''[x]+intertVlocal[x] y3[x]+x interVnl[x,B,Rb]interprod[Rb,B]==(-0.398-0.5/27.211) y3[x],y3[0.01]==f3z[0.01],y3'[0.01]==f3z'[0.01]},y3,{x,0.01,10}];

root3=FindRoot[y3[10] /. sol3[B],{B,interminBRb[Rb]}];

valBRbm[[i]]={Rb,B/.root3};

,{i,maxi}]

Show[ListContourPlot[countour2,PlotRange->{0,5},ColorFunction->"Rainbow"],ListPlot[{Rb/.fit,B/.fit},PlotStyle->Green],ListPlot[valBRbM,PlotStyle->{Orange,PointSize[0.01]}],ListPlot[valBRbm,PlotStyle->{Yellow,PointSize[0.01]}],ListPlot[valBRb,PlotStyle->{Red}],Plot[fB[Rb],{Rb,minRb,maxRb},PlotStyle->Yellow],FrameLabel->{"Rb","B"}]