Examine meson wavfunctions.
Spin parallel/antiparallel case

\(We\ define\ \(q\_parallel\) \((x)\)\ to\ be\ the\ probability\ of\ finding\ a\ quark\ with\ momentum\)
fraction x and the same spin as the antiquark.

\(We\ define\ \(q\_anti\) \((x)\)\ to\ be\ the\ probability\ of\ finding\ a\ quark\ with\ momentum\)
\(fraction\ x\ and\ the\ opposite\ spin\ as\ the\ \(\(antiquark\)\(.\)\)\)

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 spin parallel/antiparallel quark 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.

RunThrough["gunzip -c meson_5_40+4_5.m.gz", ""];
kk = Apply[Plus, Map[Last, basis[[1]]]];
pp = Table[{i/kk, 0}, {i, 1/2, kk - 1/2}];
aa = Table[{i/kk, 0}, {i, 1/2, kk - 1/2}]; 
    Do[If[basis[[i, 1, 1]] == 
          basis[[i, -1, 1]], 
        pp[[basis[[i, 1, -1]] + 1/2, 2]] += 
          Abs[wave[0][[i]]]^2 basis[[i, 
                1, -1]]/2; 
        pp[[basis[[i, -1, -1]] + 1/2, 2]] += 
          Abs[wave[
                    0][[i]]]^2 basis[[i, -1, -1]]/2, 
        aa[[basis[[i, 1, -1]] + 1/2, 2]] += 
          Abs[wave[0][[i]]]^2 basis[[i, 
                1, -1]]/2; 
        aa[[basis[[i, -1, -1]] + 1/2, 2]] += 
          Abs[wave[
                    0][[i]]]^2 basis[[i, -1, -1]]/2], {i, Length[basis]}]; 
data[kk, parallel] = pp; 
data[kk, anti] = aa;

Output to data file.

Definition[data] >> "meson_5+4_5_structure.m"

This is the proper normalization:

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

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

Compute spin parallel/antiparallel 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]]]];
pp = Table[0, {i, 1/2, kk - 1/2}, {i, 1/2, kk - 1/2}];
aa = Table[0, {i, 1/2, kk - 1/2}, {i, 1/2, kk - 1/2}]; 
    Do[If[basis[[i, 1, 1]] == basis[[i, -1, 1]], 
        pp[[basis[[i, 1, -1]] + 1/2, 
              basis[[i, -1, -1]] + 1/2]] += 
          Abs[wave[0][[i]]]^2 kk^2/2; 
        pp[[basis[[i, -1, -1]] + 1/2, 
              basis[[i, 1, -1]] + 1/2]] += 
          Abs[wave[0][[i]]]^2 kk^2/2, 
        aa[[basis[[i, 1, -1]] + 1/2, 
              basis[[i, -1, -1]] + 1/2]] += 
          Abs[wave[0][[i]]]^2 kk^2/2; 
        aa[[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, parallel] = pp;
twod[kk, anti] = aa;

Output to data file.

Definition[data] >> "meson_5+4_5_twod.m"

This is the normalization (the following should add to 1):

totalpara = Apply[Plus, pp, {0, 1}]/kk^2

totalanti = Apply[Plus, aa, {0, 1}]/kk^2

totalpara + totalanti

This is not so useful:

Block[{range = {{1/(2 kk), 1 - 1/(2 kk)}, {1/(2 kk), 
              1 - 1/(2 kk)}}}, 
      Show[ListPlot3D[pp, MeshRange -> range, ClipFill -> None, 
          DisplayFunction -> Identity], 
        ListPlot3D[aa, MeshRange -> range, ClipFill -> None, 
          DisplayFunction -> Identity], 
        ViewPoint -> {1.0, -1.4, 2.6}, 
        DisplayFunction -> $DisplayFunction, 
        PlotRange -> {0.1, 2}]]

Find the fractional spin parallel/antiparallel asymmetry

percent = 
        Array[If[
              aa[[##]] > 
                0, (aa[[##]] - pp[[##]])/(aa[[##]] + 
                    pp[[##]]), 10] &, Dimensions[aa]];

pdist = 
      Block[{range = {{1/(2 kk), 1 - 1/(2 kk)}, {1/(2 kk), 
                1 - 1/(2 kk)}}}, 
        Show[ListPlot3D[percent, MeshRange -> range, 
            ClipFill -> None, 
            AxesLabel -> {"x", "y", ""}, 
            TextStyle -> {FontSize -> 14, 
                FontFamily -> "Times"}, 
            DisplayFunction -> Identity], 
          ViewPoint -> {0.9, -1.4, 0.6}, PlotRange -> {0, 0.6}, 
          DisplayFunction -> $DisplayFunction]]

Export["relative_20.eps", pdist, ImageSize -> 6*72]

Analyze spin parallel/antiparallel quark distributions

<< "meson_5+4_5_structure.m";

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

str = data[10, parallel];

fitting = 
      NonLinearFit[Drop[Drop[str, 1], -1], 
        x^beta*(a + b*Sqrt[x] + c x)*(1 - x)^alpha, 
        x, {{alpha, 0.19}, {beta, 0.001}, {a, 0.1}, {b, -0.12}, {c, 
            0.4}}]

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

Another arbitraray choice:

fitting = 
      NonLinearFit[
        Drop[Drop[str, 
            2], -1], (a +  b LegendreP[1, 2 x - 1] +  
              c LegendreP[2, 2 x - 1]) x^e (1 - x)^d, 
        x, {{a, 0.12}, {b, 0.1}, {d, 0.2}, {c, 0}, {e, 0.2}}]

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

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

Compare a fit with the actual data:

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

Do fits for all K at once (this takes some time):

Clear[sfit]; 
    Map[(sfit[#, parallel] = 
            NonLinearFit[Drop[Drop[data[#, parallel], 1], -1], 
              x^beta*(a + b*Sqrt[x] + c x)*(1 - x)^alpha, 
              x, {{alpha, 0.19}, {beta, 0.1}, {a, 0.14}, {b, -0.3}, {c, 
                  0.6}}]) &, {10, 12, 15, 20}]; 
    Map[(sfit[#, anti] = 
            NonLinearFit[Drop[Drop[data[#, anti], 1], -1], 
              x^beta*(a + b*Sqrt[x] + c x)*(1 - x)^alpha, 
              x, {{alpha, 0.19}, {beta, 0.1}, {a, 0.14}, {b, -0.3}, {c, 
                  0.6}}]) &, {10, 12, 15, 20}]

Append all the fits to the data file:

Definition[sfit] >>> "meson_5+4_5_structure.m"

Analyze spin parallel/antiparallel quark distributions

<< "meson_5+4_5_structure.m";

Extrapolate result to infinite K and divide by x:

extrapolatedparallel = 
        Table[{x, (Fit[
                    Map[{#, sfit[#, parallel]} &, {10, 12, 15, 20}], {1, 
                      1/kvalue, 1/kvalue^2}, kvalue] /. 
                  kvalue -> Infinity)/x}, {x, 0.1, 0.9999, 
            0.01}];

fitextrapolatedparallel = 
      NonLinearFit[extrapolatedparallel, 
        x^beta*(a + b*Sqrt[x] + c x)*(1 - x)^alpha, 
        x, {{alpha, 0.19}, {beta, 0.1}, {a, 0.14}, {b, -0.3}, {c, 
            0.6}}]

extrapolatedanti = 
        Table[{x, (Fit[
                    Map[{#, sfit[#, anti]} &, {10, 12, 15, 20}], {1, 
                      1/kvalue, 1/kvalue^2}, kvalue] /. 
                  kvalue -> Infinity)/x}, {x, 0.1, 0.9999, 
            0.01}];

fitextrapolatedanti = 
      NonLinearFit[extrapolatedanti, 
        x^beta*(a + b*Sqrt[x] + c x)*(1 - x)^alpha, 
        x, {{alpha, 0.19}, {beta, 0.1}, {a, 0.14}, {b, -0.3}, {c, 
            0.6}}]

Compare the fits with the K->infinity extrapolation.

Show[ListPlot[extrapolatedparallel, DisplayFunction -> Identity], 
      Plot[fitextrapolatedparallel, {x, 0.1, 1}, 
        DisplayFunction -> Identity], 
      ListPlot[extrapolatedanti, DisplayFunction -> Identity], 
      Plot[fitextrapolatedanti, {x, 0.1, 1}, 
        DisplayFunction -> Identity], 
      DisplayFunction -> $DisplayFunction]

p1 = 
      Apply[Show[
            Plot[fitextrapolatedparallel, {x, 0.1, 1}, 
              DisplayFunction -> Identity], ##, Axes -> False,  
            Frame -> True, RotateLabel -> False, 
            TextStyle -> {FontSize -> 14, 
                FontFamily -> "Times"}, PlotRange -> {0, 0.5}, 
            FrameLabel -> {"x", "q(x)"}, 
            DisplayFunction -> $DisplayFunction] &, 
        Map[{ListPlot[
                Map[{#[[1]], #[[2]]/#[[1]]} &, 
                  data[#, parallel]], 
                DisplayFunction -> Identity]} &, {10, 12, 15, 
            20}]]

p2 = 
      Apply[Show[
            Plot[fitextrapolatedanti, {x, 0.1, 1}, 
              DisplayFunction -> Identity], ##, Axes -> False,  
            Frame -> True, RotateLabel -> False, 
            TextStyle -> {FontSize -> 14, 
                FontFamily -> "Times"}, PlotRange -> {0, 0.5}, 
            FrameLabel -> {"x", "q(x)"}, 
            DisplayFunction -> $DisplayFunction] &, 
        Map[{ListPlot[
                Map[{#[[1]], #[[2]]/#[[1]]} &, data[#, anti]], 
                DisplayFunction -> Identity]} &, {10, 12, 15, 
            20}]]

Put all the fits together:

p3 = Show[p1, p2]

Export["spin_struct_parallel.eps", p1, 
      ImageSize -> 6*72]

Export["spin_struct_anti.eps", p2, 
      ImageSize -> 6*72]

Export["spin_struct_both.eps", p3, 
      ImageSize -> 6*72]