Examine meson wavfunctions.

rmserror[data_, fit_, var_] := 
      Sqrt[Apply[Plus, 
            Map[((fit /. var -> #[[1]]) - #[[2]])^2 &, 
              data]]/Length[data]]

NonLinearFit[data_, fit_, var_, coeff_List] := 
      Block[{v = If[Length[coeff] > 1, Map[First, coeff], var]}, 
        fit /. Apply[
              FindMinimum[rmserror[data, fit, var], ##, 
                  MaxIterations -> 100, 
                  AccuracyGoal -> Infinity] &, 
              coeff][[2]]]

Compute (unpolarized) quark distributions

Clear[data]

Use data files from web page.  Since they are pretty large, we don't try to save them in Mathematica.

\(\(RunThrough["gunzip -c ~/public_html/four/meson_5_16-7_1.m.gz", ""];\)\)
kk = Apply[Plus, Map[Last, basis[[1]]]]

Here, we calculate x q(x).

str = Table[{i/kk, 0}, {i, 1/2, kk - 1/2}]; 
    Do[str[[basis[[i, 1, -1]] + 1/2, 2]] += 
        Abs[wave[0][[i]]]^2 basis[[i, 1, -1]]/2; 
      str[[basis[[i, -1, -1]] + 1/2, 2]] += 
        Abs[wave[0][[i]]]^2 basis[[i, -1, -1]]/2, {i, 
        Length[basis]}];
data[kk] = str

Output to data file.

Definition[
        data] >> "~/public_html/four/meson_5-7_1_structure.m"

This is the proper normalization:

Apply[Plus, Map[(#[[-1]]/#[[1]]) &, str]]/kk

Analyze (unpolarized) quark distributions

<< "~/public_html/four/meson_5-7_1_structure.m";

str = data[21];

This is the one that Simon used (he says that it is not well-justified):

fitting = 
      NonLinearFit[Drop[Drop[str, 1], -1], 
        a*(x^b)*(1 + c*Sqrt[x])*(1 - x)^d, 
        x, {{a, 0.19}, {b, 0.001}, {c, 1.9}, {d, 0.5}}]

rmserror[Drop[Drop[str, 1], -1], fitting, x]

This works tons better:

fitting = 
      NonLinearFit[
        Drop[Drop[str, 
            1], -1], (a +  b LegendreP[1, 2 x - 1] +  
              c LegendreP[2, 2 x - 1]) (1 - x)^d, 
        x, {{a, 0.14}, {b, 2.77}, {c, 1}, {d, 0.5}}]

rmserror[Drop[Drop[str, 1], -1], fitting, x]

Do all fits at once:

Clear[sfit]; 
    Map[(sfit[#] = 
            NonLinearFit[
              Drop[Drop[data[#], 
                  1], -1], (a +  b LegendreP[1, 2 x - 1] +  
                    c LegendreP[2, 2 x - 1]) (1 - x)^d, 
              x, {{a, 0.14}, {b, 2.77}, {c, 1}, {d, 0.5}}]) &, {8, 11, 15, 
        21}]

Append all the fits to the data file:

Definition[
        sfit] >>> "~/public_html/four/meson_5-7_1_structure.m"

Show[ListPlot[str, DisplayFunction -> Identity], 
      Plot[fitting, {x, 0, 1}, DisplayFunction -> Identity], 
      Axes -> False,  Frame -> True, 
      DisplayFunction -> $DisplayFunction]

Extrapolate result to infinite K:

Plot[
      Fit[Map[{#, sfit[#]} &, {8, 11, 15, 21}], {1, 1/kvalue, 1/kvalue^2}, 
          kvalue] /. kvalue -> Infinity, {x, 0, 1}]

Put all the fits together:

Apply[
      Show[Plot[
            Fit[Map[{#, sfit[#]} &, {8, 11, 15, 21}], {1, 1/kvalue, 
                  1/kvalue^2}, kvalue] /. kvalue -> Infinity, {x, 0, 1}, 
            DisplayFunction -> Identity], ##, Axes -> False,  
          Frame -> True, RotateLabel -> False, 
          TextStyle -> {FontSize -> 14, 
              FontFamily -> "Times"}, 
          FrameLabel -> {"x", "x q(x)"}, 
          DisplayFunction -> $DisplayFunction] &, 
      Map[{ListPlot[data[#], DisplayFunction -> Identity], 
            Plot[sfit[#], {x, 0, 1}, 
              DisplayFunction -> Identity]} &, {8, 11, 15, 
          21}]]

rho ((Subscript[J, z])=1 component)unpolarized distribution function x q(x)
\(Three\ link\ truncation, \ 
    K = 8, 11, 15, 
    21, 
along\ with\ extrapolation\ to\ K = \(\(∞\)\(.\)\)\)

Pion distribution function x q(x)
\(Three\ link\ truncation, \ 
    K = 8, 11, 15, 
    21, 
along\ with\ extrapolation\ to\ K = \(\(∞\)\(.\)\)\)

Width 600, left at zero to print

Apply[
      Show[Plot[(Fit[
                    Map[{#, sfit[#]} &, {8, 11, 15, 21}], {1, 1/kvalue, 
                      1/kvalue^2}, kvalue] /. kvalue -> Infinity)/
              x, {x, 0, 1}, DisplayFunction -> Identity], ##, 
          Axes -> False,  Frame -> True, 
          PlotRange -> {0, 1.5}, RotateLabel -> False, 
          TextStyle -> {FontSize -> 14, 
              FontFamily -> "Times"}, 
          FrameLabel -> {"x", "q(x)"}, 
          DisplayFunction -> $DisplayFunction] &, 
      Map[{ListPlot[Map[{#[[1]], #[[2]]/#[[1]]} &, data[#]], 
              DisplayFunction -> Identity], 
            Plot[sfit[#]/x, {x, 0, 1}, 
              DisplayFunction -> Identity]} &, {8, 11, 15, 
          21}]]

Pion distribution function q(x)
\(Three\ link\ truncation, \ 
    K = 8, 11, 15, 
    21, 
along\ with\ extrapolation\ to\ K = \(\(∞\)\(.\)\)\)

Width 600, left at zero to print

Compute polarized quark distributions

Clear[datap, datam]

Use data files from web page.  Since they are pretty large, we don't try to save them in Mathematica.

\(\(RunThrough["gunzip -c ~/public_html/four/meson_5_16-7_1.m.gz", ""];\)\)
kk = Apply[Plus, Map[Last, basis[[1]]]]

Calculate x q(X)

strp = Table[{i/kk, 0}, {i, 1/2, kk - 1/2}]; 
    strm = Table[{i/kk, 0}, {i, 1/2, kk - 1/2}]; 
    Do[If[basis[[i, 1, 1]] == 1, 
        strp[[basis[[i, 1, -1]] + 1/2, 2]] += 
          Abs[wave[0][[i]]]^2 basis[[i, 1, -1]]/2, 
        strm[[basis[[i, 1, -1]] + 1/2, 2]] += 
          Abs[wave[0][[i]]]^2 basis[[i, 1, -1]]/2]; 
      If[basis[[i, -1, 1]] == 1, 
        strp[[basis[[i, -1, -1]] + 1/2, 2]] += 
          Abs[wave[0][[i]]]^2 basis[[i, -1, -1]]/2, 
        strm[[basis[[i, -1, -1]] + 1/2, 2]] += 
          Abs[wave[0][[i]]]^2 basis[[i, -1, -1]]/2], {i, 
        Length[basis]}];
datap[kk] = strp;  datam[kk] = strm;

This is the proper normalization:

(Apply[Plus, Map[(#[[-1]]/#[[1]]) &, strp]] + 
          Apply[Plus, 
            Map[(#[[-1]]/#[[1]]) &, strm]])/
      kk

Definition[datap, 
        datam] >> "~/public_html/four/meson_5-7_1_polarized.m"

Analyze polarized quark distributions for (Subscript[J, z])=1

<< "~/public_html/four/meson_5-7_1_polarized.m";

Do all fits at once:

Clear[sfitp, 
      sfitm]; {Map[(sfitp[#] = 
              NonLinearFit[
                Drop[Drop[datap[#], 
                    1], -1], (a +  b LegendreP[1, 2 x - 1] +  
                      c LegendreP[2, 2 x - 1] +  
                      e LegendreP[3, 2 x - 1]) (1 - x)^d, 
                x, {{a, 0.1}, {b, -0.05}, {c, 0.1}, {d, 0.5}, {e, 
                    0.01}}]) &, {8, 11, 15, 21}], 
      Map[(sfitm[#] = 
              NonLinearFit[
                Drop[Drop[datam[#], 
                    1], -1], (a +  b LegendreP[1, 2 x - 1] +  
                      c LegendreP[2, 2 x - 1] +  
                      e LegendreP[3, 2 x - 1]) (1 - x)^d, 
                x, {{a, 0.2}, {b, 0.1}, {c, -0.1}, {d, 0.5}, {e, 
                    0.01}}]) &, {8, 11, 15, 21}]}

Append all the fits to the data file:

Definition[sfitp, 
        sfitm] >>> "~/public_html/four/meson_5-7_1_polarized.m"

Put all the fits together:

Apply[
      Show[Plot[
            Fit[Map[{#, sfitp[#]} &, {8, 11, 15, 21}], {1, 1/kvalue, 
                  1/kvalue^2}, kvalue] /. kvalue -> Infinity, {x, 0, 1}, 
            DisplayFunction -> Identity], 
          Plot[Fit[
                Map[{#, sfitm[#]} &, {8, 11, 15, 21}], {1, 1/kvalue, 
                  1/kvalue^2}, kvalue] /. kvalue -> Infinity, {x, 0, 1}, 
            DisplayFunction -> Identity], ##, Axes -> False,  
          Frame -> True, RotateLabel -> False, 
          TextStyle -> {FontSize -> 14, 
              FontFamily -> "Times"}, 
          FrameLabel -> {"x", "x q(x)"}, 
          DisplayFunction -> $DisplayFunction] &, 
      Map[{ListPlot[datap[#], DisplayFunction -> Identity, 
              PlotStyle -> {Hue[0]}], 
            ListPlot[datam[#], DisplayFunction -> Identity, 
              PlotStyle -> {Hue[2/3]}], 
            Plot[sfitp[#], {x, 0, 1}, DisplayFunction -> Identity, 
              PlotStyle -> {Hue[0]}], 
            Plot[sfitm[#], {x, 0, 1}, DisplayFunction -> Identity, 
              PlotStyle -> {Hue[2/3]}]} &, {8, 11, 15, 21}]]

rho ((Subscript[J, z])=1 component)polarized distribution functions x q(x)
\(Three\ link\ truncation, \ 
    K = 8, 11, 15, 
    21, 
along\ with\ extrapolation\ to\ K = \(\(∞\)\(.\)\)\)

Width 600, left at zero to print

Compute 2-d distributions

Clear[data, parallel, anti]

Use data files from web page.  Since they are pretty large, we don't try to save them in Mathematica.
The idea here is to cycle through the various values of K (by hand) and then dump the results for all K in a data file.

datadir = "~/public_html/four/";

RunThrough["gunzip -c " <> 
          datadir <> "meson_5_40+4_5.m.gz", ""];
kk = Apply[Plus, Map[Last, basis[[1]]]];
dist = Table[0, {i, 1/2, kk - 1/2}, {i, 1/2, kk - 1/2}];
Do[
        dist[[basis[[i, 1, -1]] + 1/2, 
              basis[[i, -1, -1]] + 1/2]] += 
          Abs[wave[0][[i]]]^2 kk^2/2; 
        dist[[basis[[i, -1, -1]] + 1/2, 
              basis[[i, 1, -1]] + 1/2]] += 
          Abs[wave[0][[i]]]^2 kk^2/2, {i, 
          Length[basis]}];
twod[kk] = dist;

This is the normalization:

total = Apply[Plus, dist, {0, 1}]/kk^2

ptwo = 
      Block[{range = {{1/(2 kk), 1 - 1/(2 kk)}, {1/(2 kk), 
                1 - 1/(2 kk)}}}, 
        Show[ListPlot3D[Map[If[# == 0, -0.1, #] &, dist, {2}], 
            MeshRange -> range, ClipFill -> None, 
            DisplayFunction -> Identity], PlotRange -> {0, 4}, 
          ViewPoint -> {2.4, -0.2, 1.8}, 
          TextStyle -> {FontSize -> 14, 
              FontFamily -> "Times"}, 
          DisplayFunction -> $DisplayFunction, 
          AxesLabel -> {"x", "y", "q(x,y)"}]]

Export["2d_struct.eps", ptwo, ImageSize -> 6*72]