data=Import["C:\\Users\\Natalia Morais\\MEOCloud\\problema10e11\\findCarbono\\pot_c_vofr.dat","Table"];

pseudoP=Table[{data[[i]][[1]],data[[i]][[3]]},{i,Dimensions[data][[1]]}];
screening=Table[{data[[i]][[1]],data[[i]][[4]]},{i,Dimensions[data][[1]]}];

local=Table[{data[[i]][[1]],pseudoP[[i]][[2]]+screening[[i]][[2]]},{i,Dimensions[data][[1]]}];

maxplot=Position[local,10.0563015`][[1]][[1]];
plotlocal=Take[local,maxplot];
ListPlot[plotlocal,PlotRange->All,AxesLabel->{"r","Vlocal"}]

interlocal=Interpolation[local]

Plot[interlocal[r],{r,0,20},PlotRange->All,AxesLabel->{"r","Vlocal"}]

localfour[q_]:= 4Pi NIntegrate[interlocal[r] r^2 SphericalBesselJ[0,q r],{r,0,20}]

Dynamic[k]
tlocalfour=Table[{k, localfour[k]},{k,1,20,0.2}];

ListPlot[tlocalfour,PlotRange->{-5,1},AxesLabel->{"k","Vlocal"}]

intertlocalfour=Interpolation[tlocalfour]

Show[Plot[intertlocalfour[q],{q,0,20},PlotRange->{-5,1}],ListPlot[tlocalfour,PlotRange->{-5,1}],AxesLabel->{"k","Vlocal"}]

Pseudo[k_,Ra_,ktf_,qz_]:=16π(((k^2)/(qz^2))-1)((E^(-(((k^2)(Ra^2))/4)))/((k^2)+(ktf^2)))

rt=FindRoot[intertlocalfour[q],{q,6}]

qfit=q/.rt;

model[k_,Ra_,ktf_]=Pseudo[k,Ra,ktf,qfit]

fit=FindFit[tlocalfour,model[k,Ra,ktf],{{Ra,0.93},{ktf,0.2}},k]

Show[Plot[model[k,Ra,ktf]/.fit,{k,0,20},PlotRange->{-5,1}],ListPlot[tlocalfour,PlotRange->{-5,1}],AxesLabel->{"k","Vlocal"}]

Plot[{intertlocalfour[q],model[q,Ra/.fit,ktf/.fit]},{q,0,20},PlotRange->{-5,1},AxesLabel->{"k","Vlocal"},PlotLegends->{"from reference","our model"}]

min[Ra_,ktf_]:=NIntegrate[(((intertlocalfour[q]-model[q,Ra,ktf]))^2)(q^2),{q,1,10}]

Dynamic[Ra]
Dynamic[ktf]
minRa=0.18;
maxRa=0.30;
dRa=0.005;
minktf=0.5;
maxktf=0.8;
dktf=0.01;
contour=Table[{Ra,ktf,min[Ra,ktf]},{Ra,minRa,maxRa,dRa},{ktf,minktf,maxktf,dktf}];

contour2=Flatten[contour,1];

Show[ListContourPlot[contour2,GridLines->Automatic],ListPlot[{{Ra/.fit,ktf/.fit}},PlotStyle->Green],FrameLabel->{"Ra","ktf"}]

siV[r_,Ra_,ktf_]:=NIntegrate[model[k,Ra,ktf]k^2 SphericalBesselJ[0,k r],{k,0,25}]/(2 Pi^2)

Dynamic[r]
Dynamic[Ra]
Dynamic[ktf]
tsiV=Table[{{r,Ra,ktf},siV[r,Ra,ktf]},{r,0,20,0.5},{Ra,minRa,maxRa,dRa},{ktf,minktf,maxktf,0.05}];

datasiV=Flatten[tsiV,2];
intersiV=Interpolation[datasiV]

valRaktf=Table[,{Ra,minRa,maxRa,dRa}];
maxi=Dimensions[valRaktf][[1]];

Dynamic[Ra]
Do[

Ra=minRa+(i-1)*dRa;

sol2[ktf_]:=NDSolve[{-(1/2)  y2''[x]+intersiV[x,Ra,ktf]  y2[x]+((1(1+1))/(2(x^2)))y2[x]==-0.199y2[x],y2[0.01]==(0.01^2),y2'[0.01]==2*(0.01^1)},y2,{x,0.01,20}];

root2=FindRoot[y2[19.9] /. sol2[ktf],{ktf,0.8}];

valRaktf[[i]]={Ra,ktf/.root2}

,{i,maxi}]
Clear[Ra]

ListPlot[valRaktf,PlotStyle->{Red,PointSize[0.01]},PlotRange->All,AxesLabel->{"Ra","ktf"}]

valRaktfb=Take[valRaktf,Dimensions[valRaktf][[1]]-1]

ListPlot[valRaktfb,PlotStyle->{Red,PointSize[0.01]},PlotRange->All,AxesLabel->{"Ra","ktf"}]

fktf[Ra_]=Fit[valRaktfb,{1,Ra,Ra^2},Ra]

Show[ListPlot[valRaktfb,PlotStyle->{Red,PointSize[0.01]},PlotRange->All],Plot[fktf[Ra],{Ra,minRa,maxRa},PlotStyle->Yellow],FrameLabel->{"Ra","ktf"}]

Show[ListContourPlot[contour2,GridLines->Automatic,ColorFunction->"Rainbow"],ListPlot[valRaktf,PlotStyle->{Red,PointSize[0.01]}],ListPlot[{{Ra/.fit,ktf/.fit}},PlotStyle->Green],Plot[fktf[Ra],{Ra,minRa,maxRa},PlotStyle->Yellow],FrameLabel->{"Ra","ktf"}]

valRaktfM=Table[,{Ra,minRa,maxRa,dRa}];
maxi=Dimensions[valRaktf][[1]];

Dynamic[Ra]
Dynamic[r]
Do[

Ra=minRa+(i-1)*dRa;

sol2[ktf_]:=NDSolve[{-(1/2)  y2''[x]+intersiV[x,Ra,ktf]  y2[x]+((1(1+1))/(2(x^2)))y2[x]==(-0.199+0.5/27.211)y2[x],y2[0.01]==(0.01^2),y2'[0.01]==2*(0.01^1)},y2,{x,0.01,20}];

root2=FindRoot[y2[10] /. sol2[ktf],{ktf,0.8}];

valRaktfM[[i]]={Ra,ktf/.root2}

,{i,maxi}]
Clear[Ra]

valRaktfm=Table[,{Ra,minRa,maxRa,dRa}];
maxi=Dimensions[valRaktf][[1]];

Dynamic[Ra]
Dynamic[r]
Do[

Ra=minRa+(i-1)*dRa;

sol2[ktf_]:=NDSolve[{-(1/2)  y2''[x]+intersiV[x,Ra,ktf]  y2[x]+((1(1+1))/(2(x^2)))y2[x]==(-0.199-0.5/27.211)y2[x],y2[0.01]==(0.01^2),y2'[0.01]==2*(0.01^1)},y2,{x,0.01,20}];

root2=FindRoot[y2[10] /. sol2[ktf],{ktf,0.8}];

valRaktfm[[i]]={Ra,ktf/.root2}

,{i,maxi}]
Clear[Ra]

Show[ListContourPlot[contour2,GridLines->Automatic,ColorFunction->"Rainbow"],ListPlot[valRaktf,PlotStyle->{Red,PointSize[0.01]}],ListPlot[valRaktfM,PlotStyle->{Orange,PointSize[0.01]}],ListPlot[valRaktfm,PlotStyle->{Yellow,PointSize[0.01]}],ListPlot[{{Ra/.fit,ktf/.fit}},PlotStyle->Green],Plot[fktf[Ra],{Ra,minRa,maxRa},PlotStyle->Yellow],FrameLabel->{"Ra","ktf"}]