(Kramer Kroenig EQ Huang, et al.,)
Hello, does anyone know how to integrate the function below to solve for the Phase Amp using Kramers Kroenig Relation's in python. The Integral can be definite through the range of 1.2 - 8 however I am unsure on how to properly integrate the singularity point using the quad function. I have tried specifying singularity points with 'points' and weights like with 'cauchys' I am unsure if I am doing it correctly as the integral goes wild. R(E') takes values from a text file with values for the points of the integral.
def reflection_function(energy):
return np.interp(energy, ref_ener, ref_ref)
def integrand(e_prime, energy):
return np.log(reflection_function(e_prime))/(e_prime**2 - energy**2)
integral_result = quad(integrand, 1, 8, args=(energy,), weight = 'cauchy', wvar = energy)[0]
print(integral_result)
It returns wild values like -220310228.0672309 and it should be around -0.7655680665661697. I am unsure why this is happening. Also I get this error
IntegrationWarning: The maximum number of subdivisions (50) has been achieved.
If increasing the limit yields no improvement it is advised to analyze
the integrand in order to determine the difficulties. If the position of a
local difficulty can be determined (singularity, discontinuity) one will
probably gain from splitting up the interval and calling the integrator
on the subranges. Perhaps a special-purpose integrator should be used.
integral_result = quad(integrand, 1, 8, args=(energy,), weight =
Putting
weight="cauchy", wvar=energy
will automatically apply a weighting1/(eprime-energy)
associated with this pole, so this doesn't need to be included in the integrand.Thus, you can keep your statement using quad, but you need to send it a modified integrand:
(This does assume that your reflection_function is positive, but since we can't see it we can't do anything about that.)
I'm also assuming that your energy is positive - use the other pole if not.
As a simpler example, the following program would find the principal value of the integral of 1/(x2-1) from 0 to 3. The integrand has a simple pole at x=1, so is modified in the call to quad by removing the factor (x-1) from the denominator. The principal value at the end should be (1/2)ln(1/2), or -0.34657...