Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Problem with a little big ray parameter #7

Open
yangjudiao opened this issue Aug 5, 2021 · 1 comment
Open

Problem with a little big ray parameter #7

yangjudiao opened this issue Aug 5, 2021 · 1 comment

Comments

@yangjudiao
Copy link

yangjudiao commented Aug 5, 2021

hello akuhara!
when I try to invert my own data setting a little bit large ray parameter like 0.15(s/km) in the param file, it just crashed and report an error like :
Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.
but when using values around 0.06-0.1(s/km) like those in the examples, it works well.
Is it a problem of the propagator matrix itself or I set something else wrong?
yangjudiao

@SuwenJunliu
Copy link

Hi Yangjudiao,

I checked the program according to your parameter, and I am pretty sure that it was caused by the e_inverse function in forward.f90

    e_inv(:,:) = (0.d0, 0.d0)
    eta = sqrt(1.d0/(beta*beta) - p*p)
    xi  = sqrt(1.d0/(alpha*alpha) - p*p)
    bp = 1.d0 - 2.d0*beta*beta*p*p

If the p is too large, the 1.d0/(beta*beta) - p*p will become a negative number.
You could check the Aki & Richards, pp. 158, Eq. (5.65), or just the hint by Akuhara Aki & Richards, pp. 161, Eq. (5.71), maybe just change all those eta and xi to complex number

    real(8) :: bp
    !real(8) :: eta, xi, bp
    complex(8) :: eta, xi
    
    e_inv(:,:) = (0.d0, 0.d0)
    eta = sqrt(cmplx(1.d0/(beta*beta) - p*p, 0))
    xi  = sqrt(cmplx(1.d0/(alpha*alpha) - p*p, 0))
    bp = 1.d0 - 2.d0*beta*beta*p*p

But do we really need to deal with those much larger rayp?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants