<<Statistics`NonlinearFit`
<<NumericalMath`SplineFit`
<<NumericalMath`Approximations`
(*________________________________________________________________KLICK____*)

(*__________INPUT (your data)_remove_comment_lines_________________KLICK____*)

(*datensatz = "/users/ferrero/webdir/det_optfun/lsf.asc";*) 

(*________READ IN DATA_(remove_comment_lines)_________________________________*)

(*tab1 = ReadList[datensatz,{Number,Number}];
peak = Max[#[[2]]& /@ tab1];
neg = Min[#[[2]]& /@ tab1];
Print[":Number of points in data set ",Length[tab1]]
plotlist = ListPlot[tab1,PlotRange->{neg,peak+peak/10},PlotStyle->{PointSize[0.02]},
AxesLabel->{"w","Exp.(MTF)"}]
mtferg = #[[2]]& /@ tab1;*)

(*______________INPUT__ test data (comment out next three instructions when you input your own data__KLICK____*)

tab1 = {{0, 0.5285133602086343}, {0.5, 0.4787424181744356},
 {1, 0.3744007289223427}, {1.5, 0.2680045857672853}, 
 {2, 0.1845856439359268}, {2.5, 0.1306828827454481}, 
 {3, 0.1046312441033665}};
plotlist = ListPlot[tab1,PlotRange->{0,0.55},PlotStyle->{PointSize[0.02]}]
mtferg = #[[2]]& /@ tab1;

integ = 15 ;    (* number of points in the PSF function  *)
q = 0.2     ;    (* upper limit of r for the PSF function *)
                    
(*_______________________SPLINE FIT________________________________________*)

halb = Length[tab1]
fourtab = tab1;                                   
spline = SplineFit[fourtab,Cubic];

(*_____________________READ OUT OF THE POLYNOMIALS__________________________*)

poly = halb-1;   
coefs = Table[i,{i,1,poly}];   
Do[coefs[[i]] = spline[[4,i]],{i,1,poly}]
polys = Table[i,{i,1,poly}];
Do[polys[[i]] = #[[1]]+t#[[2]]+t^2#[[3]]+t^3#[[4]]& /@coefs[[i]],{i,1,poly}]
ywert = Table[i,{i,1,poly}];
Do[zw = polys[[i]];
          ywert[[i]] =zw[[2]],{i,1,poly}]
          
(*splineplot = ParametricPlot[spline[u],{u,0,halb-1},AxesLabel->{"w","Spline fit(MTF)"},
PlotRange->All,Compiled->False,PlotStyle->{RGBColor[0,1,0.]}]
Show [plotlist,splineplot]*)

(*_______________HANKEL TRANSFORM __(MTF->PSF)_(might take a long time)_____*)

psferg = Table[i,{i,1,integ}];
rvect = Table[i-1,{i,1,q+1,q/integ}];
Do[psferg[[j]] = 0 ,{j,1,integ}]
Do[
Do[psferg[[j]] = (2.*Pi)*NIntegrate[ywert[[i]] BesselJ[0,2*Pi(t+i-1)*
rvect[[j]]]*(t+i-1),
{t,0,1}]+psferg[[j]]
,{i,1,poly}]
,{j,1,integ}]
psferg;
res = Table[{rvect[[i]],psferg[[i]]},{i,1,integ}];
plotpsf = ListPlot[res,PlotRange->{-2,psferg[[1]]+psferg[[1]]/10},AxesLabel->{"r","PSF"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]},PlotJoined->True]

(*____________________________________SCALING ____________________________*)

scal =(tab1[[2,1]]-tab1[[1,1]]);
tabpsf = Table[{rvect[[i]]/scal,psferg[[i]]*scal^2},{i,1,integ}];
plotpsfgen = ListPlot[tabpsf,PlotRange->{-2,(psferg[[1]]+psferg[[1]]/10)*scal^2},
AxesLabel->{"r","PSF scaled"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]},PlotJoined->True]

(*_______FOURIER TRANSFORM__MTF->LSF_____________________KLICK_____________*)

four = InverseFourier[mtferg];
Length[four];
foumax = Max[Abs[four]];
fourpkt = Table[{i-1,Abs[four[[i]]]},{i,1,Length[four]}];                              
plot = ListPlot[Abs[fourpkt],PlotRange->{0,foumax},
AxesLabel->{"x","Abs(LSF)"},  \
PlotStyle->{PointSize[0.05],RGBColor[1,0.,0.]},PlotJoined->True]
(*plot = ListPlot[Abs[fourpkt],PlotRange->{0,foumax},
AxesLabel->{"x","Abs(LSF)"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]}]*)

(*______________________________SCALING__________________________________*)

dohm = scal;
(* df = scal; *)
df = dohm/(2.0 Pi);
n = Length[mtferg];
dx = 1/(n*df);
wvect = Table[i,{i,0,(Length[fourpkt]-1)*dx,dx}];
tablsf = Table[{wvect[[i]],Abs[fourpkt[[i,2]]]},{i,1,Length[fourpkt]}];
plot = ListPlot[tablsf,AxesLabel->{"x","Abs(LSF)"},  
PlotStyle->{PointSize[0.03],RGBColor[1,0.,0.]},PlotJoined->True] ;

(*____________SAVE DATA____________________________KLICK_______*)

writeMatrixToSpreadSheet[filename_String,data_List] :=
  Module[{file=OpenWrite[filename]},
  Scan[
       (
        WriteString[file,FortranForm[First[#]]];
        Scan[
             WriteString[file,"\t",FortranForm[#]]&,
             Rest[#]
        ];
        WriteString[file,"\n"]
       )&,
        data
      ];
      Close[file]
  ]     

tabpsf;
writeMatrixToSpreadSheet["myfile.dat",tabpsf]
tablsf;
writeMatrixToSpreadSheet["myfile2.dat",tablsf]