（个人学习笔记，仅供参考）

import numpy as np
from scipy.special import kv, erfc
from scipy.integrate import dblquad
import matplotlib.pyplot as plt
import scipy.integrate as spi

w = 0.6198
g0_sq = 21.5989
rho = 0.782908
Q2 = 17.8
P2 = 0.01
z0 = 6.0

def compute_integral(x):result, error = dblquad(integrand, 1.0e-6, 50.0, lambda z: 1.0e-6, lambda z: 13.0, args=(x, g0_sq, rho, Q2, P2), epsabs=1.0e-8, epsrel=1.0e-8)print(f"Integral result: {result}, Error: {error}")  # Debug outputreturn result, error# Define range of x values (use a logarithmic scale)
x_values = np.logspace(-5, -1, 100)  # Logarithmic range for x# Initialize arrays to store results
integral_results178 = []
errors = []# Loop over x_values to compute the integral for each x
for x in x_values:result, error = compute_integral(x)integral_results178.append(result)errors.append(error)# Convert results to numpy arrays for further processing or plotting
integral_results178 = np.array(integral_results178)
errors = np.array(errors)

x178 = [0.0033, 0.0166] 
x178_err = [0, 0]
result178 = [0.428, 0.295]
result178_err = [0.061, 0.019] 

plt.figure(figsize=(5, 3))
plt.plot(x_values, integral_results178, label='$Q^2=17.8$GeV$^2$', color='blue')
plt.errorbar(x178, result178, yerr=result178_err, xerr=x178_err, fmt='d', color='black', capsize=6, label='OPAL.Data')
plt.xlabel('x')
plt.xscale('log')
plt.xlim(1e-5, 0.1)
plt.ylim(0, 2)  # Adjust y-axis based on max result
plt.ylabel('$F_2^{\gamma}(x,Q^2)$')
plt.legend()
plt.show()