how can i write a code to solve kramers-kronig equation?

hi everybody. I want to use kramers-kronig relation. I know the value of k in the wavelength range of 300-2500 nm (omega in the kramers-kronig relation can be calculated using wavelength). In this relation P means principal value of the integral. can anybody help me to write a code for this equation? please find the equation in attached figure I have tried to write a code but it does not give the correct answer

 clear all
close all
clc
A=pi;
Data = xlsread('Data.xlsx', 'Sheet1');
w = Data(:,1);
v= Data(:,1);
k= Data(:,2);
for j=1:800
     r=0;
    for i=1:800
        if i==j

        else
           y=(k(i)*(w(i+1)-w(i))/(w(i)-v(j)));
           r=r+y;
        end
    end
    z(j)=r/A;
    n(j)=1+2*z(j);
end
xlswrite('refractive index.xlsx',n',1);

[ramezan sahebi]

figure

Answer 1

[Rémi Capillon]

There are several potential sources of error here, but my first question would be, how wrong is the answer? My second one would be about the experimental data and what it looks like. I think you should be able to obtain decent results with the code you've presented. I've successfully calculated viscoelastic moduli using a similar method, but a lot of potential problems should be considered. All in all, the numerical integration method, aside from obvious bugs I could have missed, would be the last thing I would check.

Possible sources of error include :

1) The frequency range of your experimental data is too narrow ;

2) The sampling of your experimental data is not good enough ;

3) The singularity in the Kramers-Kronig relation isn't treated well enough.

I'll try to give some more details about these problems.

About 1). You'll notice in your equation that you should integrate k(omega) over the frequency range [0,+inf[, which means that a lot of information on the behaviour of k(omega) may be missing from your experimental wavelength range of [300,2500] nm. Most notably, we can see that you lack the knowledge of low frequency behaviour. Lower frequencies may be playing a vital role ! Equivalently, higher frequencies may play an equivalently important role.

About 2). The sampling of the experimental data is also of crucial importance. Your data may be smooth, but the equation does not get that information put in if the spacing of your experimental data is not small enough or not regular enough.

About 3). Function f(omega,omega') = k(omega')/(omega' - omega) explodes at the singularity and its neighbourhood. While you can't evaluate the function for omega = omega', it is of importance to give the integral some information about the behaviour of f(omega,omega') close to the singularity. I don't know if this is a problem or not because you didn't give any information about the sampling.

Now, what I would recommend to address each of these issues.

For 1), this really depends on your experimental measurement capabilities. If you cannot change the wavelength range of measurement then just leave it out and hope for the best. If you can do something about it, I would recommend including low frequency behaviour as much as possible. Widening the frequency range's upper bound could help also, but it's probably a constraint. There are few articles dealing with problems like this but a potential solution would be to use knowledge of the asymptotic behaviour of k(omega) and add it to your experimental data. The frequency (or wavelenght) spacing may be different for this part of the data and not as good, and may require more complex integration techniques for which I will give details if you are interested, but it should not be your primary concern right now. Fix the rest first.

For 2), if k(omega) is smooth enough, try to interpolate the experimental data over a much larger number of frequency points with a constant and as small as possible step.

For 3), this is also closely related to the sampling of the experimental data. The smaller step the better since it will bring you that much closer to the singularity point in each of the integrals you perform. Try to plot the function I called f(omega,omega') for a few values of omega, ideally at least a small and a large one (as small and as large as you can provide). This will help figuring out if the singularity is taken into account well enough in your numerical integration.

Good luck to you and please give any additional details you may have on the data and methodology!

Answer 2

[Shaily_T]

There is a helpful toolbox and function which help you calculate the real and imaginary parts of refractive index by K-K relations.

This is the link to it: https://www.mathworks.com/matlabcentral/fileexchange/8135-tools-for-data-analysis-in-optics-acoustics-signal-processing

Also, this book is helpful and it has some useful codes at the end: Kramers-Kronig Relations in Optical Materials Research

Hope this helps!