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

(*___________INPUT____________________________________________KLICK______*)


datensatz = "psf.asc";


(*____________READ IN DATA__________________________________________*)

tab1 = ReadList[datensatz,{Number,Number}];
peak = Max[#[[2]]& /@ tab1];
neg = Min[#[[2]]& /@ tab1];
Print["Number of points : ",Length[tab1]]
plotlist = ListPlot[tab1,PlotRange->{neg,peak+peak/10}]
qwer = Position[#[[2]]& /@ tab1,Max[#[[2]]& /@ tab1]]-1;
qwerty = qwer[[1]];
punktre =qwerty[[1]];

(*___________READ OUT OFFSET, MAXIMUM POSITION_________________________*)

qwer = Position[#[[2]]& /@ tab1,Max[#[[2]]& /@ tab1]]-1;
qwerty = qwer[[1]];
punktre =qwerty[[1]];
offset = tab1[[punktre+1,1]];
Print["Maximum after ",punktre," points"]
Print["Offset : ",offset]
dx = Abs[tab1[[1,1]]-tab1[[2,1]]];
Print["dx = ",dx]

(*____________________INPUT_____________________________KLICK___________*)



punktli = 236; (*no of points to be dropped out on the right side of the peak *)




(*_________Eliminate non necessary points_______________________________*)

tab1li= Drop[tab1,-punktli];
tab1re = Drop[tab1li,punktre];

(*__________________OFFSET,SCALING____________________________________*)

tab1b = Table[i,{i,1,Length[tab1re]}];
tab2 = Table[i,{i,1,Length[tab1re]}];
Do[tab1b[[i]] = tab1re[[i]]-{offset,0},{i,1,Length[tab1re]}];
plotlistunsk = ListPlot[tab1b,PlotRange->{neg,peak+peak/10},
PlotStyle->{PointSize[0.05],RGBColor[1,0.,0.]},
AxesLabel->{"r","PSF Exp.,without Offset"}]
scal = 1/(tab1b[[2,1]]-tab1b[[1,1]]);
Print["Number of points : ",Length[tab1b]]
Print["Scaling factor: ",scal]
Do[tab2[[i]] = tab1b[[i]]*{scal,1},{i,1,Length[tab1re]}];
tab2a=#[[2]]& /@ tab2;
n = Length[tab2a];
plotlist = ListPlot[tab2,PlotRange->{neg,peak+peak/10},
PlotStyle->{PointSize[0.05],RGBColor[1,0.,0.]},
AxesLabel->{"r","PSF,scaled"}]

(*______________INPUT____________________________KLICK_________________*)

            
integ = 10  ;            (* Number of points in MTF *)
q = 5       ;            (* upper limit of r for MTF  *)
                    
(*_______________________SPLINE FIT______________________________________*)

halb = n   ;
rvect = Table[i,{i,1,halb}];     
spline = SplineFit[tab2,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}]

(*_________HANKEL TRANSFORM_____(PSF->MTF)_________________*)

mtferg = Table[i,{i,1,integ}];
rvect = Table[i-1,{i,1,q+1,q/integ}];
Do[mtferg[[j]] = 0 ,{j,1,integ}]
Do[
Do[mtferg[[j]] = (2.*Pi)^(-1)*NIntegrate[ywert[[i]] BesselJ[0,(t+i-1)*
rvect[[j]]]*(t+i-1),
{t,0,1}]+mtferg[[j]]
,{i,1,poly}]
,{j,1,integ}]
res = Table[{rvect[[i]],mtferg[[i]]},{i,1,integ}];
(*splineplot = ParametricPlot[spline[u],{u,0,halb-1},
AxesLabel->{"r","Splinefit(PSF,scaled)"},
PlotRange->All,Compiled->False,PlotStyle->{RGBColor[0,1,0.]}]
Show[plotlist,splineplot]*)
plotmtf = ListPlot[res,PlotRange->{0,mtferg[[1]]+mtferg[[1]]/10},
AxesLabel->{"w","MTF"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]}]

(*_________________________________SCALING_______________________________*)

tabmtf = Table[{rvect[[i]]/dx,mtferg[[i]]*dx^2},{i,1,integ}];
plotmtf = ListPlot[tabmtf,PlotRange->{0,mtferg[[1]]*dx^2+mtferg[[1]]*dx^2/10},
AxesLabel->{"w","MTF"},  \
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),PKT"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]}]

(*_________________________SCALING______________________________________*)

df = Abs[tabmtf[[1,1]]-tabmtf[[2,1]]];
dxlsf = 1/(integ*df);
xvect = Table[i,{i,0,(Length[tabmtf]-1)*dxlsf,dxlsf}];
tablsf = Table[{xvect[[i]],Abs[four[[i]]]},{i,1,Length[tabmtf]}];
plot = ListPlot[tablsf,PlotRange->{0,foumax},
AxesLabel->{"x","Abs(LSF),PKT"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]}]
Show[plotlistunsk]
Show[plotmtf]

(*____________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]
  ]  
tabmtf;     
writeMatrixToSpreadSheet["myfile.dat",tabmtf]

(* tabmtfmatr = MatrixForm[tabmtf]; *)
(* mtfstr = OpenWrite["mtf_daten",FormatType->OutputForm];*)
(* Write[mtfstr,tabmtfmatr]; *)
(* Close[mtfstr]             *)

tablsf;
writeMatrixToSpreadSheet["myfile2.dat",tablsf]
(*tablsfmatr = MatrixForm[tablsf];  *)
(* lsfstr = OpenWrite["lsf_daten",FormatType->OutputForm]; *)
(* Write[lsfstr,tablsfmatr]; *)
(* Close[lsfstr]         *)