(* This program transforms LSF into MTF and subsequently MTF into PSF. In order to execute the different parts of the program, klick on the lines where "KLICK" appears and type <SHIFT>+<ENTER>. Enter your own data below every "INPUT" comment line. The format of the input data set is 2 ASCII  number columns *)

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

(*__________INPUT____________________________________________KLICK_______*)


datensatz = "/users/ferrero/webdir/det_optfun/lsf.asc"; (* Name of the data set with data to be transformed *)

(*________READ IN DATA______________________________________________*)

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->{"x","Exp.(LSF)"}]

(*_________INPUT_______________________________________________KLICK______*)



punktli = 220;          (*no of points to be dropped out on the left side of the peak *)
punktre = 241;          (*no of points to be dropped out on the right side of the peak *)
offset  = 9;            (*Offset to be subtracted             *)


(*_________________Eliminate non necessary points____________________________________*)

tab1li = Drop[tab1,-punktre];
tab1re = Drop[tab1li,punktli];
tab2 = Table[i,{i,1,Length[tab1re]}];
Do[tab2[[i]] = tab1re[[i]]-{offset,0},{i,1,Length[tab1re]}];
tab2a=#[[2]]& /@ tab2;
Print["Number of points in data set: ",Length[tab2]]
Print["This number must be odd"]
plotlist = ListPlot[tab2,PlotRange->{neg,peak+peak/10},
PlotStyle->{PointSize[0.05],RGBColor[1,0.,0.]},
AxesLabel->{"x","LSF "}]
tablog = Table[{tab2[[i,1]],Log[tab2[[i,2]]]},
{i,1,Length[tab2]}];
plotlog = ListPlot[tablog,PlotStyle->{PointSize[0.05],RGBColor[1,0.,0.]},
AxesLabel->{"x","Log(LSF)"}]
dx = Abs[tab2[[1,1]]-tab2[[2,1]]];
n = Length[tab2a];

(*____________FOURIER TRANSFORM__LSF->MTF__________________________KLICK________*)


fourpktganz = Abs[Fourier[tab2a]];
z = Max[Abs[fourpkt]];
a = Min[Abs[fourpkt]];
fourpkt = Drop[fourpktganz,-(n/2-0.5)];               
df = 1/(n*dx);
dohm = 2.0 Pi df;                                           
Length[fourpkt];
(* wvect = Table[i,{i,0,(Length[fourpkt]-1)*df,df}]; *)
wvect = Table[i,{i,0,(Length[fourpkt]-1)*dohm,dohm}];
wvectspl = Table[i,{i,0,Length[fourpkt]-1}];
tabmtf = Table[{wvect[[i]],Abs[fourpkt[[i]]]},{i,1,Length[fourpkt]}];
tabmtfspl = Table[{wvectspl[[i]],Abs[fourpkt[[i]]]},{i,1,Length[fourpkt]}];
plot = ListPlot[tabmtf,AxesLabel->{"w","Abs(MTF)"},  \
PlotStyle->{PointSize[0.05],RGBColor[1,0.,0.]},PlotJoined->True]
plot = ListPlot[tabmtf,
AxesLabel->{"w","Abs(MTF),Points"},  \
PlotStyle->{PointSize[0.02],RGBColor[0.,0.,1.]}]

(*_________INPUT________________________________________KLICK_____________*)



integ = 20;             (* number of points in the PSF function  *)
q = 1.5;                (* upper limit of r for the PSF function *)
                    
(*_______________________SPLINEFIT__________________________________________*)

halb = Length[tabmtf];   
rvect = Table[i,{i,1,halb}];      
spline = SplineFit[tabmtfspl,Cubic];

(*splineplot = ParametricPlot[spline[u],{u,0,halb-1},
AxesLabel->{"w","Splinefit(MTF)"},
PlotRange->All,Compiled->False,PlotStyle->{RGBColor[0,1,0.]}]*)

(*_______________________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}]
ywert;
(* scal = df; *)
scal = dohm;

(*________________________HANKEL TRANSFORM_____(MTF->PSF)________________________*)

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}];
(*splineplot = ParametricPlot[spline[u],{u,0,halb-1},
AxesLabel->{"w","Splinefit(MTF)"},
PlotRange->All,Compiled->False,PlotStyle->{RGBColor[0,1,0.]}]*)
(*Show[plot,splineplot]*)
plotpsf = ListPlot[res,PlotRange->{-2,psferg[[1]]+psferg[[1]]/10},
AxesLabel->{"r","PSF unscaled"},  \
PlotStyle->{PointSize[0.02],RGBColor[1,0.,0.]},PlotJoined->True]

(*____________________________________SCALING____________________________*)

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]
Show[plotlist]
Show[plot]

(*____________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]   *)


tabpsf;

     
writeMatrixToSpreadSheet["myfile2.dat",tabpsf]

(* tabpsfmatr = MatrixForm[tabpsf]  *)
(* psfstr = OpenWrite["psf_daten",FormatType->OutputForm] *)
(* Write[psfstr,tabpsfmatr] *)
(* Close[psfstr]  *)