Subroutine pystfu(kf, x, q2, xpq, jbt)
  Common /hparnt/hipr1(100), ihpr2(50), hint1(100), ihnt2(50)
  Save /hparnt/
  Common /hjcrdn/yp(3, 300), yt(3, 300)
  Save /hjcrdn/
  Common /ludat1/mstu(200), paru(200), mstj(200), parj(200)
  Save /ludat1/
  Common /ludat2/kchg(500, 3), pmas(500, 4), parf(2000), vckm(4, 4)
  Save /ludat2/
  Common /pypars/mstp(200), parp(200), msti(200), pari(200)
  Save /pypars/
  Common /pyint1/mint(400), vint(400)
  Save /pyint1/
  Dimension xpq(-6:6), xq(6), tx(6), tt(6), ts(6), nehlq(8, 2), cehlq(6, 6, 2, 8, 2), cdo(3, 6, 5, 2), cow(3, 5, 4, 2)
  Data nehlq/3, 4, 7, 5, 7, 7, 7, 7, 3, 4, 7, 6, 7, 7, 7, 7/
  Data (((cehlq(ix,it,nx,1,1),ix=1,6),it=1,6), nx=1, 2)/7.677E-01, -2.087E-01, -3.303E-01, -2.517E-02, -1.570E-02, -1.000E-04, -5.326E-01, -2.661E-01, 3.201E-01, 1.192E-01, 2.434E-02, 7.620E-03, 2.162E-01, 1.881E-01, -8.375E-02, -6.515E-02, -1.743E-02, -5.040E-03, -9.211E-02, -9.952E-02, 1.373E-02, 2.506E-02, 8.770E-03, 2.550E-03, 3.670E-02, 4.409E-02, 9.600E-04, -7.960E-03, -3.420E-03, -1.050E-03, -1.549E-02, -2.026E-02, -3.060E-03, 2.220E-03, 1.240E-03, 4.100E-04, 2.395E-01, 2.905E-01, 9.778E-02, 2.149E-02, 3.440E-03, 5.000E-04, 1.751E-02, -6.090E-03, -2.687E-02, -1.916E-02, -7.970E-03, -2.750E-03, -5.760E-03, -5.040E-03, 1.080E-03, 2.490E-03, 1.530E-03, 7.500E-04, 1.740E-03, 1.960E-03, 3.000E-04, -3.400E-04, -2.900E-04, -1.800E-04, -5.300E-04, -6.400E-04, -1.700E-04, 4.000E-05, 6.000E-05, 4.000E-05, 1.700E-04, 2.200E-04, 8.000E-05, 1.000E-05, -1.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,1,2),ix=1,6),it=1,6), nx=1, 2)/7.237E-01, -2.189E-01, -2.995E-01, -1.909E-02, -1.477E-02, 2.500E-04, -5.314E-01, -2.425E-01, 3.283E-01, 1.119E-01, 2.223E-02, 7.070E-03, 2.289E-01, 1.890E-01, -9.859E-02, -6.900E-02, -1.747E-02, -5.080E-03, -1.041E-01, -1.084E-01, 2.108E-02, 2.975E-02, 9.830E-03, 2.830E-03, 4.394E-02, 5.116E-02, -1.410E-03, -1.055E-02, -4.230E-03, -1.270E-03, -1.991E-02, -2.539E-02, -2.780E-03, 3.430E-03, 1.720E-03, 5.500E-04, 2.410E-01, 2.884E-01, 9.369E-02, 1.900E-02, 2.530E-03, 2.400E-04, 1.765E-02, -9.220E-03, -3.037E-02, -2.085E-02, -8.440E-03, -2.810E-03, -6.450E-03, -5.260E-03, 1.720E-03, 3.110E-03, 1.830E-03, 8.700E-04, 2.120E-03, 2.320E-03, 2.600E-04, -4.900E-04, -3.900E-04, -2.300E-04, -6.900E-04, -8.200E-04, -2.000E-04, 7.000E-05, 9.000E-05, 6.000E-05, 2.400E-04, 3.100E-04, 1.100E-04, 0.000E+00, -2.000E-05, -2.000E-05/
  Data (((cehlq(ix,it,nx,2,1),ix=1,6),it=1,6), nx=1, 2)/3.813E-01, -8.090E-02, -1.634E-01, -2.185E-02, -8.430E-03, -6.200E-04, -2.948E-01, -1.435E-01, 1.665E-01, 6.638E-02, 1.473E-02, 4.080E-03, 1.252E-01, 1.042E-01, -4.722E-02, -3.683E-02, -1.038E-02, -2.860E-03, -5.478E-02, -5.678E-02, 8.900E-03, 1.484E-02, 5.340E-03, 1.520E-03, 2.220E-02, 2.567E-02, -3.000E-05, -4.970E-03, -2.160E-03, -6.500E-04, -9.530E-03, -1.204E-02, -1.510E-03, 1.510E-03, 8.300E-04, 2.700E-04, 1.261E-01, 1.354E-01, 3.958E-02, 8.240E-03, 1.660E-03, 4.500E-04, 3.890E-03, -1.159E-02, -1.625E-02, -9.610E-03, -3.710E-03, -1.260E-03, -1.910E-03, -5.600E-04, 1.590E-03, 1.590E-03, 8.400E-04, 3.900E-04, 6.400E-04, 4.900E-04, -1.500E-04, -2.900E-04, -1.800E-04, -1.000E-04, -2.000E-04, -1.900E-04, 0.000E+00, 6.000E-05, 4.000E-05, 3.000E-05, 7.000E-05, 8.000E-05, 2.000E-05, -1.000E-05, -1.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,2,2),ix=1,6),it=1,6), nx=1, 2)/3.578E-01, -8.622E-02, -1.480E-01, -1.840E-02, -7.820E-03, -4.500E-04, -2.925E-01, -1.304E-01, 1.696E-01, 6.243E-02, 1.353E-02, 3.750E-03, 1.318E-01, 1.041E-01, -5.486E-02, -3.872E-02, -1.038E-02, -2.850E-03, -6.162E-02, -6.143E-02, 1.303E-02, 1.740E-02, 5.940E-03, 1.670E-03, 2.643E-02, 2.957E-02, -1.490E-03, -6.450E-03, -2.630E-03, -7.700E-04, -1.218E-02, -1.497E-02, -1.260E-03, 2.240E-03, 1.120E-03, 3.500E-04, 1.263E-01, 1.334E-01, 3.732E-02, 7.070E-03, 1.260E-03, 3.400E-04, 3.660E-03, -1.357E-02, -1.795E-02, -1.031E-02, -3.880E-03, -1.280E-03, -2.100E-03, -3.600E-04, 2.050E-03, 1.920E-03, 9.800E-04, 4.400E-04, 7.700E-04, 5.400E-04, -2.400E-04, -3.900E-04, -2.400E-04, -1.300E-04, -2.600E-04, -2.300E-04, 2.000E-05, 9.000E-05, 6.000E-05, 4.000E-05, 9.000E-05, 1.000E-04, 2.000E-05, -2.000E-05, -2.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,3,1),ix=1,6),it=1,6), nx=1, 2)/6.870E-02, -6.861E-02, 2.973E-02, -5.400E-03, 3.780E-03, -9.700E-04, -1.802E-02, 1.400E-04, 6.490E-03, -8.540E-03, 1.220E-03, -1.750E-03, -4.650E-03, 1.480E-03, -5.930E-03, 6.000E-04, -1.030E-03, -8.000E-05, 6.440E-03, 2.570E-03, 2.830E-03, 1.150E-03, 7.100E-04, 3.300E-04, -3.930E-03, -2.540E-03, -1.160E-03, -7.700E-04, -3.600E-04, -1.900E-04, 2.340E-03, 1.930E-03, 5.300E-04, 3.700E-04, 1.600E-04, 9.000E-05, 1.014E+00, -1.106E+00, 3.374E-01, -7.444E-02, 8.850E-03, -8.700E-04, 9.233E-01, -1.285E+00, 4.475E-01, -9.786E-02, 1.419E-02, -1.120E-03, 4.888E-02, -1.271E-01, 8.606E-02, -2.608E-02, 4.780E-03, -6.000E-04, -2.691E-02, 4.887E-02, -1.771E-02, 1.620E-03, 2.500E-04, -6.000E-05, 7.040E-03, -1.113E-02, 1.590E-03, 7.000E-04, -2.000E-04, 0.000E+00, -1.710E-03, 2.290E-03, 3.800E-04, -3.500E-04, 4.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,3,2),ix=1,6),it=1,6), nx=1, 2)/1.008E-01, -7.100E-02, 1.973E-02, -5.710E-03, 2.930E-03, -9.900E-04, -5.271E-02, -1.823E-02, 1.792E-02, -6.580E-03, 1.750E-03, -1.550E-03, 1.220E-02, 1.763E-02, -8.690E-03, -8.800E-04, -1.160E-03, -2.100E-04, -1.190E-03, -7.180E-03, 2.360E-03, 1.890E-03, 7.700E-04, 4.100E-04, -9.100E-04, 2.040E-03, -3.100E-04, -1.050E-03, -4.000E-04, -2.400E-04, 1.190E-03, -1.700E-04, -2.000E-04, 4.200E-04, 1.700E-04, 1.000E-04, 1.081E+00, -1.189E+00, 3.868E-01, -8.617E-02, 1.115E-02, -1.180E-03, 9.917E-01, -1.396E+00, 4.998E-01, -1.159E-01, 1.674E-02, -1.720E-03, 5.099E-02, -1.338E-01, 9.173E-02, -2.885E-02, 5.890E-03, -6.500E-04, -3.178E-02, 5.703E-02, -2.070E-02, 2.440E-03, 1.100E-04, -9.000E-05, 8.970E-03, -1.392E-02, 2.050E-03, 6.500E-04, -2.300E-04, 2.000E-05, -2.340E-03, 3.010E-03, 5.000E-04, -3.900E-04, 6.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,4,1),ix=1,6),it=1,6), nx=1, 2)/9.482E-01, -9.578E-01, 1.009E-01, -1.051E-01, 3.456E-02, -3.054E-02, -9.627E-01, 5.379E-01, 3.368E-01, -9.525E-02, 1.488E-02, -2.051E-02, 4.300E-01, -8.306E-02, -3.372E-01, 4.902E-02, -9.160E-03, 1.041E-02, -1.925E-01, -1.790E-02, 2.183E-01, 7.490E-03, 4.140E-03, -1.860E-03, 8.183E-02, 1.926E-02, -1.072E-01, -1.944E-02, -2.770E-03, -5.200E-04, -3.884E-02, -1.234E-02, 5.410E-02, 1.879E-02, 3.350E-03, 1.040E-03, 2.948E+01, -3.902E+01, 1.464E+01, -3.335E+00, 5.054E-01, -5.915E-02, 2.559E+01, -3.955E+01, 1.661E+01, -4.299E+00, 6.904E-01, -8.243E-02, -1.663E+00, 1.176E+00, 1.118E+00, -7.099E-01, 1.948E-01, -2.404E-02, -2.168E-01, 8.170E-01, -7.169E-01, 1.851E-01, -1.924E-02, -3.250E-03, 2.088E-01, -4.355E-01, 2.239E-01, -2.446E-02, -3.620E-03, 1.910E-03, -9.097E-02, 1.601E-01, -5.681E-02, -2.500E-03, 2.580E-03, -4.700E-04/
  Data (((cehlq(ix,it,nx,4,2),ix=1,6),it=1,6), nx=1, 2)/2.367E+00, 4.453E-01, 3.660E-01, 9.467E-02, 1.341E-01, 1.661E-02, -3.170E+00, -1.795E+00, 3.313E-02, -2.874E-01, -9.827E-02, -7.119E-02, 1.823E+00, 1.457E+00, -2.465E-01, 3.739E-02, 6.090E-03, 1.814E-02, -1.033E+00, -9.827E-01, 2.136E-01, 1.169E-01, 5.001E-02, 1.684E-02, 5.133E-01, 5.259E-01, -1.173E-01, -1.139E-01, -4.988E-02, -2.021E-02, -2.881E-01, -3.145E-01, 5.667E-02, 9.161E-02, 4.568E-02, 1.951E-02, 3.036E+01, -4.062E+01, 1.578E+01, -3.699E+00, 6.020E-01, -7.031E-02, 2.700E+01, -4.167E+01, 1.770E+01, -4.804E+00, 7.862E-01, -1.060E-01, -1.909E+00, 1.357E+00, 1.127E+00, -7.181E-01, 2.232E-01, -2.481E-02, -2.488E-01, 9.781E-01, -8.127E-01, 2.094E-01, -2.997E-02, -4.710E-03, 2.506E-01, -5.427E-01, 2.672E-01, -3.103E-02, -1.800E-03, 2.870E-03, -1.128E-01, 2.087E-01, -6.972E-02, -2.480E-03, 2.630E-03, -8.400E-04/
  Data (((cehlq(ix,it,nx,5,1),ix=1,6),it=1,6), nx=1, 2)/4.968E-02, -4.173E-02, 2.102E-02, -3.270E-03, 3.240E-03, -6.700E-04, -6.150E-03, -1.294E-02, 6.740E-03, -6.890E-03, 9.000E-04, -1.510E-03, -8.580E-03, 5.050E-03, -4.900E-03, -1.600E-04, -9.400E-04, -1.500E-04, 7.840E-03, 1.510E-03, 2.220E-03, 1.400E-03, 7.000E-04, 3.500E-04, -4.410E-03, -2.220E-03, -8.900E-04, -8.500E-04, -3.600E-04, -2.000E-04, 2.520E-03, 1.840E-03, 4.100E-04, 3.900E-04, 1.600E-04, 9.000E-05, 9.235E-01, -1.085E+00, 3.464E-01, -7.210E-02, 9.140E-03, -9.100E-04, 9.315E-01, -1.274E+00, 4.512E-01, -9.775E-02, 1.380E-02, -1.310E-03, 4.739E-02, -1.296E-01, 8.482E-02, -2.642E-02, 4.760E-03, -5.700E-04, -2.653E-02, 4.953E-02, -1.735E-02, 1.750E-03, 2.800E-04, -6.000E-05, 6.940E-03, -1.132E-02, 1.480E-03, 6.500E-04, -2.100E-04, 0.000E+00, -1.680E-03, 2.340E-03, 4.200E-04, -3.400E-04, 5.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,5,2),ix=1,6),it=1,6), nx=1, 2)/6.478E-02, -4.537E-02, 1.643E-02, -3.490E-03, 2.710E-03, -6.700E-04, -2.223E-02, -2.126E-02, 1.247E-02, -6.290E-03, 1.120E-03, -1.440E-03, -1.340E-03, 1.362E-02, -6.130E-03, -7.900E-04, -9.000E-04, -2.000E-04, 5.080E-03, -3.610E-03, 1.700E-03, 1.830E-03, 6.800E-04, 4.000E-04, -3.580E-03, 6.000E-05, -2.600E-04, -1.050E-03, -3.800E-04, -2.300E-04, 2.420E-03, 9.300E-04, -1.000E-04, 4.500E-04, 1.700E-04, 1.100E-04, 9.868E-01, -1.171E+00, 3.940E-01, -8.459E-02, 1.124E-02, -1.250E-03, 1.001E+00, -1.383E+00, 5.044E-01, -1.152E-01, 1.658E-02, -1.830E-03, 4.928E-02, -1.368E-01, 9.021E-02, -2.935E-02, 5.800E-03, -6.600E-04, -3.133E-02, 5.785E-02, -2.023E-02, 2.630E-03, 1.600E-04, -8.000E-05, 8.840E-03, -1.416E-02, 1.900E-03, 5.800E-04, -2.500E-04, 1.000E-05, -2.300E-03, 3.080E-03, 5.500E-04, -3.700E-04, 7.000E-05, 1.000E-05/
  Data (((cehlq(ix,it,nx,6,1),ix=1,6),it=1,6), nx=1, 2)/9.270E-03, -1.817E-02, 9.590E-03, -6.390E-03, 1.690E-03, -1.540E-03, 5.710E-03, -1.188E-02, 6.090E-03, -4.650E-03, 1.240E-03, -1.310E-03, -3.960E-03, 7.100E-03, -3.590E-03, 1.840E-03, -3.900E-04, 3.400E-04, 1.120E-03, -1.960E-03, 1.120E-03, -4.800E-04, 1.000E-04, -4.000E-05, 4.000E-05, -3.000E-05, -1.800E-04, 9.000E-05, -5.000E-05, -2.000E-05, -4.200E-04, 7.300E-04, -1.600E-04, 5.000E-05, 5.000E-05, 5.000E-05, 8.098E-01, -1.042E+00, 3.398E-01, -6.824E-02, 8.760E-03, -9.000E-04, 8.961E-01, -1.217E+00, 4.339E-01, -9.287E-02, 1.304E-02, -1.290E-03, 3.058E-02, -1.040E-01, 7.604E-02, -2.415E-02, 4.600E-03, -5.000E-04, -2.451E-02, 4.432E-02, -1.651E-02, 1.430E-03, 1.200E-04, -1.000E-04, 1.122E-02, -1.457E-02, 2.680E-03, 5.800E-04, -1.200E-04, 3.000E-05, -7.730E-03, 7.330E-03, -7.600E-04, -2.400E-04, 1.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,6,2),ix=1,6),it=1,6), nx=1, 2)/9.980E-03, -1.945E-02, 1.055E-02, -6.870E-03, 1.860E-03, -1.560E-03, 5.700E-03, -1.203E-02, 6.250E-03, -4.860E-03, 1.310E-03, -1.370E-03, -4.490E-03, 7.990E-03, -4.170E-03, 2.050E-03, -4.400E-04, 3.300E-04, 1.470E-03, -2.480E-03, 1.460E-03, -5.700E-04, 1.200E-04, -1.000E-05, -9.000E-05, 1.500E-04, -3.200E-04, 1.200E-04, -6.000E-05, -4.000E-05, -4.200E-04, 7.600E-04, -1.400E-04, 4.000E-05, 7.000E-05, 5.000E-05, 8.698E-01, -1.131E+00, 3.836E-01, -8.111E-02, 1.048E-02, -1.300E-03, 9.626E-01, -1.321E+00, 4.854E-01, -1.091E-01, 1.583E-02, -1.700E-03, 3.057E-02, -1.088E-01, 8.022E-02, -2.676E-02, 5.590E-03, -5.600E-04, -2.845E-02, 5.164E-02, -1.918E-02, 2.210E-03, -4.000E-05, -1.500E-04, 1.311E-02, -1.751E-02, 3.310E-03, 5.100E-04, -1.200E-04, 5.000E-05, -8.590E-03, 8.380E-03, -9.200E-04, -2.600E-04, 1.000E-05, -1.000E-05/
  Data (((cehlq(ix,it,nx,7,1),ix=1,6),it=1,6), nx=1, 2)/9.010E-03, -1.401E-02, 7.150E-03, -4.130E-03, 1.260E-03, -1.040E-03, 6.280E-03, -9.320E-03, 4.780E-03, -2.890E-03, 9.100E-04, -8.200E-04, -2.930E-03, 4.090E-03, -1.890E-03, 7.600E-04, -2.300E-04, 1.400E-04, 3.900E-04, -1.200E-03, 4.400E-04, -2.500E-04, 2.000E-05, -2.000E-05, 2.600E-04, 1.400E-04, -8.000E-05, 1.000E-04, 1.000E-05, 1.000E-05, -2.600E-04, 3.200E-04, 1.000E-05, -1.000E-05, 1.000E-05, -1.000E-05, 8.029E-01, -1.075E+00, 3.792E-01, -7.843E-02, 1.007E-02, -1.090E-03, 7.903E-01, -1.099E+00, 4.153E-01, -9.301E-02, 1.317E-02, -1.410E-03, -1.704E-02, -1.130E-02, 2.882E-02, -1.341E-02, 3.040E-03, -3.600E-04, -7.200E-04, 7.230E-03, -5.160E-03, 1.080E-03, -5.000E-05, -4.000E-05, 3.050E-03, -4.610E-03, 1.660E-03, -1.300E-04, -1.000E-05, 1.000E-05, -4.360E-03, 5.230E-03, -1.610E-03, 2.000E-04, -2.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,7,2),ix=1,6),it=1,6), nx=1, 2)/8.980E-03, -1.459E-02, 7.510E-03, -4.410E-03, 1.310E-03, -1.070E-03, 5.970E-03, -9.440E-03, 4.800E-03, -3.020E-03, 9.100E-04, -8.500E-04, -3.050E-03, 4.440E-03, -2.100E-03, 8.500E-04, -2.400E-04, 1.400E-04, 5.300E-04, -1.300E-03, 5.600E-04, -2.700E-04, 3.000E-05, -2.000E-05, 2.000E-04, 1.400E-04, -1.100E-04, 1.000E-04, 0.000E+00, 0.000E+00, -2.600E-04, 3.200E-04, 0.000E+00, -3.000E-05, 1.000E-05, -1.000E-05, 8.672E-01, -1.174E+00, 4.265E-01, -9.252E-02, 1.244E-02, -1.460E-03, 8.500E-01, -1.194E+00, 4.630E-01, -1.083E-01, 1.614E-02, -1.830E-03, -2.241E-02, -5.630E-03, 2.815E-02, -1.425E-02, 3.520E-03, -4.300E-04, -7.300E-04, 8.030E-03, -5.780E-03, 1.380E-03, -1.300E-04, -4.000E-05, 3.460E-03, -5.380E-03, 1.960E-03, -2.100E-04, 1.000E-05, 1.000E-05, -4.850E-03, 5.950E-03, -1.890E-03, 2.600E-04, -3.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,8,1),ix=1,6),it=1,6), nx=1, 2)/4.410E-03, -7.480E-03, 3.770E-03, -2.580E-03, 7.300E-04, -7.100E-04, 3.840E-03, -6.050E-03, 3.030E-03, -2.030E-03, 5.800E-04, -5.900E-04, -8.800E-04, 1.660E-03, -7.500E-04, 4.700E-04, -1.000E-04, 1.000E-04, -8.000E-05, -1.500E-04, 1.200E-04, -9.000E-05, 3.000E-05, 0.000E+00, 1.300E-04, -2.200E-04, -2.000E-05, -2.000E-05, -2.000E-05, -2.000E-05, -7.000E-05, 1.900E-04, -4.000E-05, 2.000E-05, 0.000E+00, 0.000E+00, 6.623E-01, -9.248E-01, 3.519E-01, -7.930E-02, 1.110E-02, -1.180E-03, 6.380E-01, -9.062E-01, 3.582E-01, -8.479E-02, 1.265E-02, -1.390E-03, -2.581E-02, 2.125E-02, 4.190E-03, -4.980E-03, 1.490E-03, -2.100E-04, 7.100E-04, 5.300E-04, -1.270E-03, 3.900E-04, -5.000E-05, -1.000E-05, 3.850E-03, -5.060E-03, 1.860E-03, -3.500E-04, 4.000E-05, 0.000E+00, -3.530E-03, 4.460E-03, -1.500E-03, 2.700E-04, -3.000E-05, 0.000E+00/
  Data (((cehlq(ix,it,nx,8,2),ix=1,6),it=1,6), nx=1, 2)/4.260E-03, -7.530E-03, 3.830E-03, -2.680E-03, 7.600E-04, -7.300E-04, 3.640E-03, -6.050E-03, 3.030E-03, -2.090E-03, 5.900E-04, -6.000E-04, -9.200E-04, 1.710E-03, -8.200E-04, 5.000E-04, -1.200E-04, 1.000E-04, -5.000E-05, -1.600E-04, 1.300E-04, -9.000E-05, 3.000E-05, 0.000E+00, 1.300E-04, -2.100E-04, -1.000E-05, -2.000E-05, -2.000E-05, -1.000E-05, -8.000E-05, 1.800E-04, -5.000E-05, 2.000E-05, 0.000E+00, 0.000E+00, 7.146E-01, -1.007E+00, 3.932E-01, -9.246E-02, 1.366E-02, -1.540E-03, 6.856E-01, -9.828E-01, 3.977E-01, -9.795E-02, 1.540E-02, -1.790E-03, -3.053E-02, 2.758E-02, 2.150E-03, -4.880E-03, 1.640E-03, -2.500E-04, 9.200E-04, 4.200E-04, -1.340E-03, 4.600E-04, -8.000E-05, -1.000E-05, 4.230E-03, -5.660E-03, 2.140E-03, -4.300E-04, 6.000E-05, 0.000E+00, -3.890E-03, 5.000E-03, -1.740E-03, 3.300E-04, -4.000E-05, 0.000E+00/
  Data ((cdo(ip,is,1,1),is=1,6), ip=1, 3)/4.190E-01, 3.460E+00, 4.400E+00, 0.000E+00, 0.000E+00, 0.000E+00, 4.000E-03, 7.240E-01, -4.860E+00, 0.000E+00, 0.000E+00, 0.000E+00, -7.000E-03, -6.600E-02, 1.330E+00, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,1,2),is=1,6), ip=1, 3)/3.740E-01, 3.330E+00, 6.030E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.400E-02, 7.530E-01, -6.220E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, -7.600E-02, 1.560E+00, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,2,1),is=1,6), ip=1, 3)/7.630E-01, 4.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, -2.370E-01, 6.270E-01, -4.210E-01, 0.000E+00, 0.000E+00, 0.000E+00, 2.600E-02, -1.900E-02, 3.300E-02, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,2,2),is=1,6), ip=1, 3)/7.610E-01, 3.830E+00, 0.000E+00, 0.000E+00, 0.000E+00, 0.000E+00, -2.320E-01, 6.270E-01, -4.180E-01, 0.000E+00, 0.000E+00, 0.000E+00, 2.300E-02, -1.900E-02, 3.600E-02, 0.000E+00, 0.000E+00, 0.000E+00/
  Data ((cdo(ip,is,3,1),is=1,6), ip=1, 3)/1.265E+00, 0.000E+00, 8.050E+00, 0.000E+00, 0.000E+00, 0.000E+00, -1.132E+00, -3.720E-01, 1.590E+00, 6.310E+00, -1.050E+01, 1.470E+01, 2.930E-01, -2.900E-02, -1.530E-01, -2.730E-01, -3.170E+00, 9.800E+00/
  Data ((cdo(ip,is,3,2),is=1,6), ip=1, 3)/1.670E+00, 0.000E+00, 9.150E+00, 0.000E+00, 0.000E+00, 0.000E+00, -1.920E+00, -2.730E-01, 5.300E-01, 1.570E+01, -1.010E+02, 2.230E+02, 5.820E-01, -1.640E-01, -7.630E-01, -2.830E+00, 4.470E+01, -1.170E+02/
  Data ((cdo(ip,is,4,1),is=1,6), ip=1, 3)/0.000E+00, -3.600E-02, 6.350E+00, 0.000E+00, 0.000E+00, 0.000E+00, 1.350E-01, -2.220E-01, 3.260E+00, -3.030E+00, 1.740E+01, -1.790E+01, -7.500E-02, -5.800E-02, -9.090E-01, 1.500E+00, -1.130E+01, 1.560E+01/
  Data ((cdo(ip,is,4,2),is=1,6), ip=1, 3)/0.000E+00, -1.200E-01, 3.510E+00, 0.000E+00, 0.000E+00, 0.000E+00, 6.700E-02, -2.330E-01, 3.660E+00, -4.740E-01, 9.500E+00, -1.660E+01, -3.100E-02, -2.300E-02, -4.530E-01, 3.580E-01, -5.430E+00, 1.550E+01/
  Data ((cdo(ip,is,5,1),is=1,6), ip=1, 3)/1.560E+00, 0.000E+00, 6.000E+00, 9.000E+00, 0.000E+00, 0.000E+00, -1.710E+00, -9.490E-01, 1.440E+00, -7.190E+00, -1.650E+01, 1.530E+01, 6.380E-01, 3.250E-01, -1.050E+00, 2.550E-01, 1.090E+01, -1.010E+01/
  Data ((cdo(ip,is,5,2),is=1,6), ip=1, 3)/8.790E-01, 0.000E+00, 4.000E+00, 9.000E+00, 0.000E+00, 0.000E+00, -9.710E-01, -1.160E+00, 1.230E+00, -5.640E+00, -7.540E+00, -5.960E-01, 4.340E-01, 4.760E-01, -2.540E-01, -8.170E-01, 5.500E+00, 1.260E-01/
  Data ((cow(ip,is,1,1),is=1,5), ip=1, 3)/4.0000E-01, 7.0000E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -6.2120E-02, 6.4780E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -7.1090E-03, 1.3350E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00/
  Data ((cow(ip,is,1,2),is=1,5), ip=1, 3)/4.0000E-01, 6.2800E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -5.9090E-02, 6.4360E-01, 0.0000E+00, 0.0000E+00, 0.0000E+00, -6.5240E-03, 1.4510E-02, 0.0000E+00, 0.0000E+00, 0.0000E+00/
  Data ((cow(ip,is,2,1),is=1,5), ip=1, 3)/8.8800E-01, 0.0000E+00, 3.1100E+00, 6.0000E+00, 0.0000E+00, -1.8020E+00, -1.5760E+00, -1.3170E-01, 2.8010E+00, -1.7280E+01, 1.8120E+00, 1.2000E+00, 5.0680E-01, -1.2160E+01, 2.0490E+01/
  Data ((cow(ip,is,2,2),is=1,5), ip=1, 3)/7.9400E-01, 0.0000E+00, 2.8900E+00, 6.0000E+00, 0.0000E+00, -9.1440E-01, -1.2370E+00, 5.9660E-01, -3.6710E+00, -8.1910E+00, 5.9660E-01, 6.5820E-01, -2.5500E-01, -2.3040E+00, 7.7580E+00/
  Data ((cow(ip,is,3,1),is=1,5), ip=1, 3)/9.0000E-01, 0.0000E+00, 5.0000E+00, 0.0000E+00, 0.0000E+00, -2.4280E-01, -2.1200E-01, 8.6730E-01, 1.2660E+00, 2.3820E+00, 1.3860E-01, 3.6710E-03, 4.7470E-02, -2.2150E+00, 3.4820E-01/
  Data ((cow(ip,is,3,2),is=1,5), ip=1, 3)/9.0000E-01, 0.0000E+00, 5.0000E+00, 0.0000E+00, 0.0000E+00, -1.4170E-01, -1.6970E-01, -2.4740E+00, -2.5340E+00, 5.6210E-01, -1.7400E-01, -9.6230E-02, 1.5750E+00, 1.3780E+00, -2.7010E-01/
  Data ((cow(ip,is,4,1),is=1,5), ip=1, 3)/0.0000E+00, -2.2120E-02, 2.8940E+00, 0.0000E+00, 0.0000E+00, 7.9280E-02, -3.7850E-01, 9.4330E+00, 5.2480E+00, 8.3880E+00, -6.1340E-02, -1.0880E-01, -1.0852E+01, -7.1870E+00, -1.1610E+01/
  Data ((cow(ip,is,4,2),is=1,5), ip=1, 3)/0.0000E+00, -8.8200E-02, 1.9240E+00, 0.0000E+00, 0.0000E+00, 6.2290E-02, -2.8920E-01, 2.4240E-01, -4.4630E+00, -8.3670E-01, -4.0990E-02, -1.0820E-01, 2.0360E+00, 5.2090E+00, -4.8400E-02/
  alam = 0.
  Do kfl = -6, 6
    xpq(kfl) = 0.
  End Do
  If (x<0. .Or. x>1.) Then
    Write (mstu(11), 1000) x
    Return
  End If
  kfa = iabs(kf)
  If (kfa/=211 .And. kfa/=2212 .And. kfa/=2112) Then
    Write (mstu(11), 1100) kf
    Return
  End If
  If (mstp(51)==0 .Or. mstp(52)>=2) Then
    kfe = kfa
    If (kfa==2112) kfe = 2212
    Call pystfe(kfe, x, q2, xpq)
    Goto 230
  End If
  If (kfa==211) Goto 200
  If (mstp(51)==1 .Or. mstp(51)==2) Then
    nset = mstp(51)
    If (nset==1) alam = 0.2
    If (nset==2) alam = 0.29
    tmin = log(5./alam**2)
    tmax = log(1E8/alam**2)
    If (mstp(52)==0) Then
      t = tmin
    Else
      t = log(q2/alam**2)
    End If
    vt = max(-1., min(1.,(2.*t-tmax-tmin)/(tmax-tmin)))
    nx = 1
    If (x<=0.1) nx = 2
    If (nx==1) vx = (2.*x-1.1)/0.9
    If (nx==2) vx = max(-1., (2.*log(x)+11.51293)/6.90776)
    cxs = 1.
    If (x<1E-4 .And. abs(parp(51)-1.)>0.01) cxs = (1E-4/x)**(parp(51)-1.)
    tx(1) = 1.
    tx(2) = vx
    tx(3) = 2.*vx**2 - 1.
    tx(4) = 4.*vx**3 - 3.*vx
    tx(5) = 8.*vx**4 - 8.*vx**2 + 1.
    tx(6) = 16.*vx**5 - 20.*vx**3 + 5.*vx
    tt(1) = 1.
    tt(2) = vt
    tt(3) = 2.*vt**2 - 1.
    tt(4) = 4.*vt**3 - 3.*vt
    tt(5) = 8.*vt**4 - 8.*vt**2 + 1.
    tt(6) = 16.*vt**5 - 20.*vt**3 + 5.*vt
    Do kfl = 1, 6
      xqsum = 0.
      Do it = 1, 6
        Do ix = 1, 6
          xqsum = xqsum + cehlq(ix, it, nx, kfl, nset)*tx(ix)*tt(it)
        End Do
      End Do
      xq(kfl) = xqsum*(1.-x)**nehlq(kfl, nset)*cxs
    End Do
    xpq(0) = xq(4)
    xpq(1) = xq(2) + xq(3)
    xpq(2) = xq(1) + xq(3)
    xpq(3) = xq(5)
    xpq(4) = xq(6)
    xpq(-1) = xq(3)
    xpq(-2) = xq(3)
    xpq(-3) = xq(5)
    xpq(-4) = xq(6)
    If (mstp(54)>=5) Then
      If (nset==1) tmin = 8.1905
      If (nset==2) tmin = 7.4474
      If (t<=tmin) Goto 140
      vt = max(-1., min(1.,(2.*t-tmax-tmin)/(tmax-tmin)))
      tt(1) = 1.
      tt(2) = vt
      tt(3) = 2.*vt**2 - 1.
      tt(4) = 4.*vt**3 - 3.*vt
      tt(5) = 8.*vt**4 - 8.*vt**2 + 1.
      tt(6) = 16.*vt**5 - 20.*vt**3 + 5.*vt
      xqsum = 0.
      Do it = 1, 6
        Do ix = 1, 6
          xqsum = xqsum + cehlq(ix, it, nx, 7, nset)*tx(ix)*tt(it)
        End Do
      End Do
      xpq(5) = xqsum*(1.-x)**nehlq(7, nset)
      xpq(-5) = xpq(5)
      140 Continue
    End If
    If (mstp(54)>=6) Then
      If (nset==1) tmin = 11.5528
      If (nset==2) tmin = 10.8097
      tmin = tmin + 2.*log(pmas(6,1)/30.)
      tmax = tmax + 2.*log(pmas(6,1)/30.)
      If (t<=tmin) Goto 160
      vt = max(-1., min(1.,(2.*t-tmax-tmin)/(tmax-tmin)))
      tt(1) = 1.
      tt(2) = vt
      tt(3) = 2.*vt**2 - 1.
      tt(4) = 4.*vt**3 - 3.*vt
      tt(5) = 8.*vt**4 - 8.*vt**2 + 1.
      tt(6) = 16.*vt**5 - 20.*vt**3 + 5.*vt
      xqsum = 0.
      Do it = 1, 6
        Do ix = 1, 6
          xqsum = xqsum + cehlq(ix, it, nx, 8, nset)*tx(ix)*tt(it)
        End Do
      End Do
      xpq(6) = xqsum*(1.-x)**nehlq(8, nset)
      xpq(-6) = xpq(6)
      160 Continue
    End If
  Else If (mstp(51)==3 .Or. mstp(51)==4) Then
    nset = mstp(51) - 2
    If (nset==1) alam = 0.2
    If (nset==2) alam = 0.4
    If (mstp(52)<=0) Then
      sd = 0.
    Else
      sd = log(log(max(q2,4.)/alam**2)/log(4./alam**2))
    End If
    Do kfl = 1, 5
      Do is = 1, 6
        ts(is) = cdo(1, is, kfl, nset) + cdo(2, is, kfl, nset)*sd + cdo(3, is, kfl, nset)*sd**2
      End Do
      If (kfl<=2) Then
        eulbt1 = pygamm(ts(1))*pygamm(ts(2)+1.)/pygamm(ts(1)+ts(2)+1.)
        xq(kfl) = x**ts(1)*(1.-x)**ts(2)*(1.+ts(3)*x)/(eulbt1*(1.+ts(3)*ts(1)/(ts(1)+ts(2)+1.)))
      Else
        xq(kfl) = ts(1)*x**ts(2)*(1.-x)**ts(3)*(1.+ts(4)*x+ts(5)*x**2+ts(6)*x**3)
      End If
    End Do
    xpq(0) = xq(5)
    xpq(1) = xq(2) + xq(3)/6.
    xpq(2) = 3.*xq(1) - xq(2) + xq(3)/6.
    xpq(3) = xq(3)/6.
    xpq(4) = xq(4)
    xpq(-1) = xq(3)/6.
    xpq(-2) = xq(3)/6.
    xpq(-3) = xq(3)/6.
    xpq(-4) = xq(4)
  Else If (mstp(51)>=11 .And. mstp(51)<=13) Then
    Call pystfe(2212, x, q2, xpq)
  Else
    Write (mstu(11), 1200) mstp(51)
  End If
  Goto 230
  200 If ((mstp(51)>=1 .And. mstp(51)<=4) .Or. (mstp(51)>=11 .And. mstp(51)<=13)) Then
    nset = 1
    If (mstp(51)==2 .Or. mstp(51)==4 .Or. mstp(51)==13) nset = 2
    If (nset==1) alam = 0.2
    If (nset==2) alam = 0.4
    If (mstp(52)<=0) Then
      sd = 0.
    Else
      sd = log(log(max(q2,4.)/alam**2)/log(4./alam**2))
    End If
    Do kfl = 1, 4
      Do is = 1, 5
        ts(is) = cow(1, is, kfl, nset) + cow(2, is, kfl, nset)*sd + cow(3, is, kfl, nset)*sd**2
      End Do
      If (kfl==1) Then
        eulbt2 = pygamm(ts(1))*pygamm(ts(2)+1.)/pygamm(ts(1)+ts(2)+1.)
        xq(kfl) = x**ts(1)*(1.-x)**ts(2)/eulbt2
      Else
        xq(kfl) = ts(1)*x**ts(2)*(1.-x)**ts(3)*(1.+ts(4)*x+ts(5)*x**2)
      End If
    End Do
    xpq(0) = xq(2)
    xpq(1) = xq(3)/6.
    xpq(2) = xq(1) + xq(3)/6.
    xpq(3) = xq(3)/6.
    xpq(4) = xq(4)
    xpq(-1) = xq(1) + xq(3)/6.
    xpq(-2) = xq(3)/6.
    xpq(-3) = xq(3)/6.
    xpq(-4) = xq(4)
  Else
    Write (mstu(11), 1200) mstp(51)
  End If
  230 If (kfa==2112) Then
    xps = xpq(1)
    xpq(1) = xpq(2)
    xpq(2) = xps
    xps = xpq(-1)
    xpq(-1) = xpq(-2)
    xpq(-2) = xps
  End If
  If (kf<0) Then
    Do kfl = 1, 4
      xps = xpq(kfl)
      xpq(kfl) = xpq(-kfl)
      xpq(-kfl) = xps
    End Do
  End If
  Do kfl = -6, 6
    xpq(kfl) = max(0., xpq(kfl))
    If (iabs(kfl)>mstp(54)) xpq(kfl) = 0.
  End Do
  If ((jbt/=1 .And. jbt/=2) .Or. ihpr2(6)==0 .Or. ihnt2(16)==1) Goto 400
  atnm = ihnt2(2*jbt-1)
  If (atnm<=1.0) Goto 400
  If (jbt==1) Then
    bbr2 = (yp(1,ihnt2(11))**2+yp(2,ihnt2(11))**2)/1.44/atnm**0.66666
  Else If (jbt==2) Then
    bbr2 = (yt(1,ihnt2(12))**2+yt(2,ihnt2(12))**2)/1.44/atnm**0.66666
  End If
  bbr2 = min(1.0, bbr2)
  abx = (atnm**0.33333333-1.0)
  apx = hipr1(6)*4.0/3.0*abx*sqrt(1.0-bbr2)
  aax = 1.192*alog(atnm)**0.1666666
  rrx = aax*(x**3-1.2*x**2+0.21*x) + 1.0 - (apx-1.079*abx*sqrt(x)/alog(atnm+1.0))*exp(-x**2.0/0.01)
  Do kfl = -6, 6
    xpq(kfl) = xpq(kfl)*rrx
  End Do
  400 Continue
  Return
  1000 Format (' Error: x value outside physical range, x =', 1P, E12.3)
  1100 Format (' Error: illegal particle code for structure function,', ' KF =', I5)
  1200 Format (' Error: bad value of parameter MSTP(51) in PYSTFU,', ' MSTP(51) =', I5)
End Subroutine pystfu
