Green's function basis

[rmodrak]

Hi all,

Thanks very much for your previous responses and for creating these very useful packages.

Recently I tried running test_get_greens_vs_get_seismograms, except rather than using a local AxiSEM database to obtain the input traces, I used the syngine web service. Strangely, the syngine synthetics do not match the source-weighted combination of syngine Green’s functions, though I copied and pasted the linear combination code directly from the instaseis test case:

I noticed that a different mapping for the moment tensor elements results in a much better agreement. The odd thing is that this new mapping (Mxx, Myy, Mzz, Mxz, Myz, Mxy = Mtt, Mpp, Mrr, Mrp, -Mrt, -Mtp) differs from the mapping at line 46 of test_greens_vs_seismograms. Here is the improved result from the new mapping (there is still a slight discrepancy on the T components, which can be fixed by making ad hoc sign changes in the linear combination):

By chance could anyone comment on the above figures? Ready-to-go scripts for generating both figures are here:

generate_figure_1.txt generate_figure_2.txt

Thanks, Ryan

[rmodrak]

Hi Lion, Thanks for your thoughts. Here are responses.

1) The scripts have been attached to the github issue.

2) Perhaps its less ambiguous to say syngine and instaseis work in the GCMT convention (up,south,east). syngine synthetics are exactly what I'd expect given this convention. There is possibly an inconsistency though involving the syngine Green's functions. I'll take a look at issue #8 .

3) Having tested locations away from the equator and poles, it's not a boundary issue.

4) I've only tried the above scripts, which use syngine queries instead of instaseis calls. (Presumably the original test_greens_vs_seismograms method works, since travis reports no errors?)

================== Hi Ryan,

I'm not sure I have time to look at this in detail until next week, but a few thoughts:

(1) You forgot to send along the attachement.

(2) Instaseis always operates in a spherical r, theta, phi system, or it should at least.

(3) The error might actually be kind of similar to your previous one as (at least in the test code) the source is again at the north pole.

(4) Does it only happen with syngine or also with local databases?

We did test this fairly extensively but there might of course be an honest bug somewhere in the line. Also please open github issues for every problem you find - otherwise it just gets lost in hour black hole like inboxes.

Cheers!

Lion

[rmodrak]

Okay, thanks to Carl Tape for just now emailing Sarah Minson.

If Sarah Minson remembers correctly, Minson & Dreger 2008 used the AkiRichards convention (north,east,down). It appears syngine uses the GCMT convention (up,south,east). The two conventions are related by

Mrr = Mzz Mtt = Mxx Mpp = Myy Mrt = Mxz Mrp = -Myz Mtp = -Mxy

Yet, even with this mapping, test_greens_vs_seismogram still doesn't work.

The remaining discrepancy can be explained by a sign error on the ZDS, RDS, and TSS Green's functions returned by syngine. Simply change the sign of these Green's functions, and everything works:

The following script was used to create the above figure: generate_figure_3.txt

I think the above comments are on the right track, but please give me a few days to double check everything...

[rmodrak]

After taking a few days to double check, I still agree with the above analysis, but I would add a few additional comments

• Sarah Minson clarified by email that (north,east,down) are the coordinate axes she had in mind in her 2008 paper. But the 2008 paper itself seems somewhat ambiguous, and it's understandable that different groups might implement the Green's functions and/or linear combination differently.

• the seiscomp definition appears inconsistent with both the (up-south-east) and (north,east,down) conventions. To my knowledge these are the two most widely used basis, being the ones employed by GCMT and Aki & Richards, respectively.

• After Carl and others double check the above comments, I'd suggest adding a defintion="AkiRichards" option for get_greens_function. This new option would differ from defintion="seiscomp" by a permutation of and sign changes on the TSS, ZSS, RSS, TDS, ZDS, RDS, ZDD, RDD, ZEP, REP traces.

[martinvandriel]

Email from Carl:

Hi Ryan, hi all,

Thanks for posting this issue/comment: https://github.com/krischer/instaseis/issues/66#issuecomment-470793928

My impression is that in order to match our MT inversion results, you need to make TWO changes: 1) assumption of north-east-down basis used in Minson and Dreger (2008) 2) sign flip on three of the green's functions (ZDS, RDS, and TSS)

I suspect that the basis issue in #1 is the source of the issue in #2, but it could be that there are two distinct issues.

To others: We can proceed with using instaseis, but we have to use these two fixes in order to get results that match our own code. I'm now "watching" instaseis so hopefully I'll see any action on this issue.

Thanks, Carl

[martinvandriel]

Hi Carl & Ryan,

sorry for beeing slow on this. I went back through my communication with Gempa when we implemented this. Our reference back then was to reproduce what they had running in production with their clients for their moment tensor tool with Green's functions from qseis.

I think both comments you make are valid:

https://github.com/krischer/instaseis/issues/8#issuecomment-111394460

https://github.com/krischer/instaseis/issues/8#issuecomment-130739969

The 'definition' keyword was readily foreseen to enable other conventions, so feel free to implement any convention you like and send a pull request. We might want to improve the documentation as well.

As with syngine, they use their own http interface and the definitions are not necessarily the same. If you query syngine via the instaseis database though db = instaseis.open_db('syngine:\\db_name'), everything should behave as a local database. If you see discrepencies there, we would need to get in touch with IRIS.

Cheers, Martin

[martinvandriel]

should have read the thread to the end. this seems to contradict Carl's point 2): https://github.com/krischer/instaseis/issues/8#issuecomment-142927825

After some discussion with @gempa-jabe decided to switch the radial sign to be consistent with Minson & Dreger. TODO!

Done in 058beb0

Can we close this issue?

[martinvandriel]

Yet another piece of confusion: https://github.com/krischer/instaseis/blob/master/instaseis/tests/test_instaseis.py#L545

    # in Minson & Dreger, 2008, z is up, so we have
    # Mtt = Mxx, Mpp = Myy, Mrr = Mzz
    # Mrp = Myz, Mrt = Mxz, Mtp = Mxy

[martinvandriel]

Now I am trying to reproduce the original problem, but somehow I don't see it: This is the code to produce the figure, adapted from the test_get_greens_vs_get_seismogram:

import matplotlib.pyplot as plt
import numpy as np
from instaseis import open_db, Source, Receiver

db = open_db('syngine://prem_a_5s')
print(db)

depth_in_m = 1000
epicentral_distance_degree = 20
azimuth = 177.

# some 'random' moment tensor
Mxx, Myy, Mzz, Mxz, Myz, Mxy = 1e20, 0.5e20, 0.3e20, 0.7e20, 0.4e20, 0.6e20

# in Minson & Dreger, 2008, z is up, so we have
# Mtt = Mxx, Mpp = Myy, Mrr = Mzz
# Mrp = Myz, Mrt = Mxz, Mtp = Mxy
source = Source(latitude=90., longitude=0., depth_in_m=depth_in_m,
                m_tt=Mxx, m_pp=Myy, m_rr=Mzz, m_rt=Mxz, m_rp=Myz, m_tp=Mxy)
receiver = Receiver(latitude=90. - epicentral_distance_degree,
                    longitude=azimuth)

st_ref = db.get_seismograms(source, receiver, components=('Z', 'R', 'T'))

st_greens = db.get_greens_function(epicentral_distance_degree,
                                   depth_in_m, definition="seiscomp")

TSS = st_greens.select(channel="TSS")[0].data  # NOQA
ZSS = st_greens.select(channel="ZSS")[0].data  # NOQA
RSS = st_greens.select(channel="RSS")[0].data  # NOQA
TDS = st_greens.select(channel="TDS")[0].data  # NOQA
ZDS = st_greens.select(channel="ZDS")[0].data  # NOQA
RDS = st_greens.select(channel="RDS")[0].data  # NOQA
ZDD = st_greens.select(channel="ZDD")[0].data  # NOQA
RDD = st_greens.select(channel="RDD")[0].data  # NOQA
ZEP = st_greens.select(channel="ZEP")[0].data  # NOQA
REP = st_greens.select(channel="REP")[0].data  # NOQA

az = np.deg2rad(azimuth)
# eq (6) in Minson & Dreger, 2008
uz = Mxx * (ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Myy * (-ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Mzz * (ZDD / 3. + ZEP / 3.) \
    + Mxy * ZSS * np.sin(2 * az) \
    + Mxz * ZDS * np.cos(az) \
    + Myz * ZDS * np.sin(az)

# eq (7) in Minson & Dreger, 2008
ur = Mxx * (RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Myy * (-RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Mzz * (RDD / 3. + REP / 3.) \
    + Mxy * RSS * np.sin(2 * az) \
    + Mxz * RDS * np.cos(az) \
    + Myz * RDS * np.sin(az)

# eq (8) in Minson & Dreger, 2008
ut = Mxx * TSS / 2. * np.sin(2 * az) \
    - Myy * TSS / 2. * np.sin(2 * az) \
    - Mxy * TSS * np.cos(2 * az) \
    + Mxz * TDS * np.sin(az) \
    - Myz * TDS * np.cos(az)

fig, ax = plt.subplots(3, 1, sharex=True)

ax[0].plot(st_ref.select(component="Z")[0].data, label='get_seismograms')
ax[0].plot(uz, label='get_greens_function')
ax[0].legend()

ax[1].plot(st_ref.select(component="R")[0].data)
ax[1].plot(ur)

ax[2].plot(st_ref.select(component="T")[0].data)
ax[2].plot(ut)

plt.show()

Is this what you where trying to do?

[rmodrak]

Let me respond to your points one at a time, in the order they were made.

I believe the above image shows the Aki & Richards convention. I think that generate_figure_3 demonstrates that defintion='seiscomp' is similar but not identical to Aki & Richards.

[rmodrak]

I agree, implementing a new convention and improving the documentation could nicely resolve everything.

[rmodrak]

Now I am trying to reproduce the original problem, but somehow I don't see it:
image
This is the code to produce the figure, adapted from the test_get_greens_vs_get_seismogram:

import matplotlib.pyplot as plt
import numpy as np
from instaseis import open_db, Source, Receiver

db = open_db('syngine://prem_a_5s')
print(db)

depth_in_m = 1000
epicentral_distance_degree = 20
azimuth = 177.

# some 'random' moment tensor
Mxx, Myy, Mzz, Mxz, Myz, Mxy = 1e20, 0.5e20, 0.3e20, 0.7e20, 0.4e20, 0.6e20

# in Minson & Dreger, 2008, z is up, so we have
# Mtt = Mxx, Mpp = Myy, Mrr = Mzz
# Mrp = Myz, Mrt = Mxz, Mtp = Mxy
source = Source(latitude=90., longitude=0., depth_in_m=depth_in_m,
                m_tt=Mxx, m_pp=Myy, m_rr=Mzz, m_rt=Mxz, m_rp=Myz, m_tp=Mxy)
receiver = Receiver(latitude=90. - epicentral_distance_degree,
                    longitude=azimuth)

st_ref = db.get_seismograms(source, receiver, components=('Z', 'R', 'T'))

st_greens = db.get_greens_function(epicentral_distance_degree,
                                   depth_in_m, definition="seiscomp")

TSS = st_greens.select(channel="TSS")[0].data  # NOQA
ZSS = st_greens.select(channel="ZSS")[0].data  # NOQA
RSS = st_greens.select(channel="RSS")[0].data  # NOQA
TDS = st_greens.select(channel="TDS")[0].data  # NOQA
ZDS = st_greens.select(channel="ZDS")[0].data  # NOQA
RDS = st_greens.select(channel="RDS")[0].data  # NOQA
ZDD = st_greens.select(channel="ZDD")[0].data  # NOQA
RDD = st_greens.select(channel="RDD")[0].data  # NOQA
ZEP = st_greens.select(channel="ZEP")[0].data  # NOQA
REP = st_greens.select(channel="REP")[0].data  # NOQA

az = np.deg2rad(azimuth)
# eq (6) in Minson & Dreger, 2008
uz = Mxx * (ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Myy * (-ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Mzz * (ZDD / 3. + ZEP / 3.) \
    + Mxy * ZSS * np.sin(2 * az) \
    + Mxz * ZDS * np.cos(az) \
    + Myz * ZDS * np.sin(az)

# eq (7) in Minson & Dreger, 2008
ur = Mxx * (RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Myy * (-RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Mzz * (RDD / 3. + REP / 3.) \
    + Mxy * RSS * np.sin(2 * az) \
    + Mxz * RDS * np.cos(az) \
    + Myz * RDS * np.sin(az)

# eq (8) in Minson & Dreger, 2008
ut = Mxx * TSS / 2. * np.sin(2 * az) \
    - Myy * TSS / 2. * np.sin(2 * az) \
    - Mxy * TSS * np.cos(2 * az) \
    + Mxz * TDS * np.sin(az) \
    - Myz * TDS * np.cos(az)

fig, ax = plt.subplots(3, 1, sharex=True)

ax[0].plot(st_ref.select(component="Z")[0].data, label='get_seismograms')
ax[0].plot(uz, label='get_greens_function')
ax[0].legend()

ax[1].plot(st_ref.select(component="R")[0].data)
ax[1].plot(ur)

ax[2].plot(st_ref.select(component="T")[0].data)
ax[2].plot(ut)

plt.show()

Is this what you where trying to do?

Yes, this is what I was trying to do in generate_figure_1 above. This generate_figure_1 script is "ready-to-go" for python2.7 in the sense that the only non-standard library dependencies are numpy, obspy, and instaseis.

The only significant differences I am seeing between your code and my generate_figure_1 are that

• your code uses a local database and mine syngine queries. However, upon further examination, I am finding that the syngine url result and open_db('syngine://name_of_db') result agree, up to an expected difference in sampling rate. So this is probably not the cause of the discrepancy.

• your code uses a source at the north pole and mine one away from the poles and equator.

I am traveling this weekend but will continue looking into the discrepancy when I get back.

[martinvandriel]

Ok, trying to move the source away from the north pole revealed the problem also in my example. It seems that everything is consistent with a single mistake: the azimuth I used in the tests was actually 180 - azimuth. That switches all signs in Minsons & Dregers formulas that contain the cos(azimuth) and sin(2 * azimuth) terms. Still need to work out how to fix that and propagate downstream.

[martinvandriel]

For documentation, here is how I think the azimuth might be the one to blame:

import matplotlib.pyplot as plt
import numpy as np
from instaseis import open_db, Source, Receiver
from obspy.geodetics.base import gps2dist_azimuth, locations2degrees

db = open_db('syngine://prem_a_5s')

depth_in_m = 1000

src_lat, src_lon = 23., 57.
rec_lat, rec_lon = 37., 17.

epicentral_distance_degree = locations2degrees(src_lat, src_lon,
                                               rec_lat, rec_lon)

azimuth = gps2dist_azimuth(src_lat, src_lon, rec_lat, rec_lon, f=0.)[1]

# FIX AZIMUTH PROBLEM
azimuth = 180. - azimuth

# some 'random' moment tensor
Mxx, Myy, Mzz, Mxz, Myz, Mxy = 1e20, 0.5e20, 0.3e20, 0.7e20, 0.4e20, 0.6e20

# assuming in Minson & Dreger, 2008, z is up, we have
# Mtt = Mxx, Mpp = Myy, Mrr = Mzz
# Mrp = Myz, Mrt = Mxz, Mtp = Mxy
source = Source(latitude=src_lat, longitude=src_lon, depth_in_m=depth_in_m,
                m_tt=Mxx, m_pp=Myy, m_rr=Mzz, m_rt=Mxz, m_rp=Myz, m_tp=Mxy)
receiver = Receiver(latitude=rec_lat, longitude=rec_lon)

st_ref = db.get_seismograms(source, receiver, components=('Z', 'R', 'T'))

st_greens = db.get_greens_function(epicentral_distance_degree,
                                   depth_in_m, definition="seiscomp")

TSS = st_greens.select(channel="TSS")[0].data
ZSS = st_greens.select(channel="ZSS")[0].data
RSS = st_greens.select(channel="RSS")[0].data
TDS = st_greens.select(channel="TDS")[0].data
ZDS = st_greens.select(channel="ZDS")[0].data
RDS = st_greens.select(channel="RDS")[0].data
ZDD = st_greens.select(channel="ZDD")[0].data
RDD = st_greens.select(channel="RDD")[0].data
ZEP = st_greens.select(channel="ZEP")[0].data
REP = st_greens.select(channel="REP")[0].data

az = np.deg2rad(azimuth)
# eq (6) in Minson & Dreger, 2008
uz = Mxx * (ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Myy * (-ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Mzz * (ZDD / 3. + ZEP / 3.) \
    + Mxy * ZSS * np.sin(2 * az) \
    + Mxz * ZDS * np.cos(az) \
    + Myz * ZDS * np.sin(az)

# eq (7) in Minson & Dreger, 2008
ur = Mxx * (RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Myy * (-RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Mzz * (RDD / 3. + REP / 3.) \
    + Mxy * RSS * np.sin(2 * az) \
    + Mxz * RDS * np.cos(az) \
    + Myz * RDS * np.sin(az)

# eq (8) in Minson & Dreger, 2008
ut = Mxx * TSS / 2. * np.sin(2 * az) \
    - Myy * TSS / 2. * np.sin(2 * az) \
    - Mxy * TSS * np.cos(2 * az) \
    + Mxz * TDS * np.sin(az) \
    - Myz * TDS * np.cos(az)

fig, ax = plt.subplots(3, 1, sharex=True)

ax[0].plot(st_ref.select(component="Z")[0].data, label='get_seismograms')
ax[0].plot(uz, label='get_greens_function')
ax[0].legend()

ax[1].plot(st_ref.select(component="R")[0].data)
ax[1].plot(ur)

ax[2].plot(st_ref.select(component="T")[0].data)
ax[2].plot(ut)

plt.show()

[martinvandriel]

And here is the verification that fixing the azimuth is consistent with the sign changes suggested by Carl and @rmodrak :

import matplotlib.pyplot as plt
import numpy as np
from instaseis import open_db, Source, Receiver
from obspy.geodetics.base import gps2dist_azimuth, locations2degrees

db = open_db('syngine://prem_a_5s')

depth_in_m = 1000

src_lat, src_lon = 23., 57.
rec_lat, rec_lon = 37., 17.

epicentral_distance_degree = locations2degrees(src_lat, src_lon,
                                               rec_lat, rec_lon)

azimuth = gps2dist_azimuth(src_lat, src_lon, rec_lat, rec_lon, f=0.)[1]

# some 'random' moment tensor
Mxx, Myy, Mzz, Mxz, Myz, Mxy = 1e20, 0.5e20, 0.3e20, 0.7e20, 0.4e20, 0.6e20

# assuming in Minson & Dreger, 2008, z is down, we have
# Mtt = Mxx, Mpp = Myy, Mrr = Mzz
# Mrp = -Myz, Mrt = Mxz, Mtp = -Mxy
source = Source(latitude=src_lat, longitude=src_lon, depth_in_m=depth_in_m,
                m_tt=Mxx, m_pp=Myy, m_rr=Mzz, m_rt=Mxz, m_rp=-Myz, m_tp=-Mxy)
receiver = Receiver(latitude=rec_lat, longitude=rec_lon)

st_ref = db.get_seismograms(source, receiver, components=('Z', 'R', 'T'))

st_greens = db.get_greens_function(epicentral_distance_degree,
                                   depth_in_m, definition="seiscomp")

TSS = -st_greens.select(channel="TSS")[0].data
ZSS = st_greens.select(channel="ZSS")[0].data
RSS = st_greens.select(channel="RSS")[0].data
TDS = st_greens.select(channel="TDS")[0].data
ZDS = -st_greens.select(channel="ZDS")[0].data
RDS = -st_greens.select(channel="RDS")[0].data
ZDD = st_greens.select(channel="ZDD")[0].data
RDD = st_greens.select(channel="RDD")[0].data
ZEP = st_greens.select(channel="ZEP")[0].data
REP = st_greens.select(channel="REP")[0].data

az = np.deg2rad(azimuth)
# eq (6) in Minson & Dreger, 2008
uz = Mxx * (ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Myy * (-ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Mzz * (ZDD / 3. + ZEP / 3.) \
    + Mxy * ZSS * np.sin(2 * az) \
    + Mxz * ZDS * np.cos(az) \
    + Myz * ZDS * np.sin(az)

# eq (7) in Minson & Dreger, 2008
ur = Mxx * (RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Myy * (-RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Mzz * (RDD / 3. + REP / 3.) \
    + Mxy * RSS * np.sin(2 * az) \
    + Mxz * RDS * np.cos(az) \
    + Myz * RDS * np.sin(az)

# eq (8) in Minson & Dreger, 2008
ut = Mxx * TSS / 2. * np.sin(2 * az) \
    - Myy * TSS / 2. * np.sin(2 * az) \
    - Mxy * TSS * np.cos(2 * az) \
    + Mxz * TDS * np.sin(az) \
    - Myz * TDS * np.cos(az)

fig, ax = plt.subplots(3, 1, sharex=True)

ax[0].plot(st_ref.select(component="Z")[0].data, label='get_seismograms')
ax[0].plot(uz, label='get_greens_function')
ax[0].legend()

ax[1].plot(st_ref.select(component="R")[0].data)
ax[1].plot(ur)

ax[2].plot(st_ref.select(component="T")[0].data)
ax[2].plot(ut)

plt.show()

[martinvandriel]

swapping ZDS, RDS and TSS is equivalent to swapping the sign of m1 and m4 here: https://github.com/krischer/instaseis/blob/master/instaseis/database_interfaces/base_instaseis_db.py#L133-L139 but I don't see yet how this fits the definitions of the coordinate systems - it would mean that the conversion from qseis to r, t,p is: r=z, t=-x, p=-y. But that is not consistent with that ascii image.

[rmodrak]

Thanks for your help getting to the bottom of this.

I agree with this comment from base_instaseis_db.py, which gives the conversion from up,south,east (r,t,p) to north,east,down (x,y,z): https://github.com/krischer/instaseis/blob/master/instaseis/database_interfaces/base_instaseis_db.py#L130-L131

If we accept that L130-L131 is the correct conversion, then this rules out the 180 - azimuth fix, because the keyword arguments given to Source differ from L130-L131.

Also, if we accept that L130-L131 is the correct conversion, then the following code should work, but it doesn't. It may be worthwhile to look carefully at L133-L152 to see if some fix is required there. I will continue looking...

import matplotlib.pyplot as plt
import numpy as np
from instaseis import open_db, Source, Receiver
from obspy.geodetics.base import gps2dist_azimuth, locations2degrees

db = open_db('syngine://prem_a_5s')

depth_in_m = 1000

src_lat, src_lon = 23., 57.
rec_lat, rec_lon = 37., 17.

epicentral_distance_degree = locations2degrees(src_lat, src_lon,
                                               rec_lat, rec_lon)

azimuth = gps2dist_azimuth(src_lat, src_lon, rec_lat, rec_lon, f=0.)[1]

# some 'random' moment tensor
Mxx, Myy, Mzz, Mxz, Myz, Mxy = 1e20, 0.5e20, 0.3e20, 0.7e20, 0.4e20, 0.6e20

# USING THE FORMULAS GIVEN BY
#https://github.com/krischer/instaseis/blob/master/instaseis/database_interfaces/base_instaseis_db.py#L130-L131
source = Source(latitude=src_lat, longitude=src_lon, depth_in_m=depth_in_m,
                m_tt=Mxx, m_pp=Myy, m_rr=Mzz, m_rt=Mxz, m_rp=-Myz, m_tp=-Mxy)
receiver = Receiver(latitude=rec_lat, longitude=rec_lon)

st_ref = db.get_seismograms(source, receiver, components=('Z', 'R', 'T'))

st_greens = db.get_greens_function(epicentral_distance_degree,
                                   depth_in_m, definition="seiscomp")

TSS = st_greens.select(channel="TSS")[0].data
ZSS = st_greens.select(channel="ZSS")[0].data
RSS = st_greens.select(channel="RSS")[0].data
TDS = st_greens.select(channel="TDS")[0].data
ZDS = st_greens.select(channel="ZDS")[0].data
RDS = st_greens.select(channel="RDS")[0].data
ZDD = st_greens.select(channel="ZDD")[0].data
RDD = st_greens.select(channel="RDD")[0].data
ZEP = st_greens.select(channel="ZEP")[0].data
REP = st_greens.select(channel="REP")[0].data

az = np.deg2rad(azimuth)
# eq (6) in Minson & Dreger, 2008
uz = Mxx * (ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Myy * (-ZSS / 2. * np.cos(2 * az) - ZDD / 6. + ZEP / 3.) \
    + Mzz * (ZDD / 3. + ZEP / 3.) \
    + Mxy * ZSS * np.sin(2 * az) \
    + Mxz * ZDS * np.cos(az) \
    + Myz * ZDS * np.sin(az)

# eq (7) in Minson & Dreger, 2008
ur = Mxx * (RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Myy * (-RSS / 2. * np.cos(2 * az) - RDD / 6. + REP / 3.) \
    + Mzz * (RDD / 3. + REP / 3.) \
    + Mxy * RSS * np.sin(2 * az) \
    + Mxz * RDS * np.cos(az) \
    + Myz * RDS * np.sin(az)

# eq (8) in Minson & Dreger, 2008
ut = Mxx * TSS / 2. * np.sin(2 * az) \
    - Myy * TSS / 2. * np.sin(2 * az) \
    - Mxy * TSS * np.cos(2 * az) \
    + Mxz * TDS * np.sin(az) \
    - Myz * TDS * np.cos(az)

fig, ax = plt.subplots(3, 1, sharex=True)

nt = 500

ax[0].plot(st_ref.select(component="Z")[0].data[:nt], label='get_seismograms')
ax[0].plot(uz[:nt], label='get_greens_function')
ax[0].legend()

ax[1].plot(st_ref.select(component="R")[0].data[:nt])
ax[1].plot(ur[:nt])

ax[2].plot(st_ref.select(component="T")[0].data[:nt])
ax[2].plot(ut[:nt])

plt.show()

[carltape]

I don't have the answer. But for me, it would help to add a couple more comments. Here's what I suggest. Even though the order of entry (Mxx, Myy, etc) may not matter in this code, it can be a source of problems in other codes.

CURRENT

# some 'random' moment tensor
Mxx, Myy, Mzz, Mxz, Myz, Mxy = 1e20, 0.5e20, 0.3e20, 0.7e20, 0.4e20, 0.6e20

# assuming in Minson & Dreger, 2008, z is down, we have
# Mtt = Mxx, Mpp = Myy, Mrr = Mzz
# Mrp = -Myz, Mrt = Mxz, Mtp = -Mxy

PROPOSED CHANGES

# some 'random' moment tensor (for the equations below)
# specified in x-y-z as north-east-down, which is assumed to be the basis in Minson & Dreger (2008)
# (and also in Aki & Richards, Jost & Herrmann 1989, Herrmann and Wang 1985)
Mxx, Myy, Mzz, Mxy, Mxz, Myz = 1e20, 0.5e20, 0.3e20, 0.6e20, 0.7e20, 0.4e20

# convert from north-east-down to up-south-east (for the Source() function)
# Mrr = Mzz, Mtt = Mxx, Mpp = Myy
# Mrt = Mxz, Mrp = -Myz, Mtp = -Mxy

The comments above match the implementation I used here: https://github.com/carltape/compearth/blob/master/momenttensor/matlab/convert_MT.m convention M = [M11 M22 M33 M12 M13 M23] Mout(1,:) = M(3,:); Mout(2,:) = M(1,:); Mout(3,:) = M(2,:); Mout(4,:) = M(5,:); Mout(5,:) = -M(6,:); Mout(6,:) = -M(4,:);

[rmodrak]

Hi Carl, It looks like your proposed changes improve the comments and make the code easier to understand, but leave the underlying mathematical formulas unchanged. So it seems Carl, myself, and L130-L131 of base_instaseis_db.py all agree regarding the mathematical formulas relating the two coordinate systems.

[martinvandriel]

then this rules out the 180 - azimuth fix

I think my wording ('fix') was misleading. What I meant was: If I use the azimuth in the same wrong way as in my original tests, then the solutions for sources anywhere are consistent with my original tests with the source at the northpole, my original understanding of the coordinate systems and your observed issues.

My original version however is obviously wrong, because it is consistent with the assumptions that (1) the sources I got from the qseis definition are in north, east, down (2) the coordinate system in Minson & Dreger is south, east, up, and (3) azimuth in Minson & Dreger eq. (6-8) is measured counterclockwise from south.

Now my missing link to the solution is this: If I assume Minson and Dreger is north, east down and azimuth measured clockwise from north (i.e. correcting assumptions (2) and (3)), than:

• either assumption (1) does not hold and the coordinates from qseis to r, t,p is r=z, t=-x, p=-y (don't believe this)
• or I got something wrong with the definitions of the sources I got from here https://github.com/krischer/instaseis/issues/8#issuecomment-111394460 (but I don't see how this could be wrong)

For me the easiest way out without breaking the seiscomp interface (assuming it's tested on their side and actually gives correct results for the moment tenors) would be to add a new definition option which only differs in the sign of m1 and m4 (https://github.com/krischer/instaseis/blob/master/instaseis/database_interfaces/base_instaseis_db.py#L133-L139). But I would be more comfortable if I understood why this is the solution. Do you have another source for the definitions of the sources assumed by Minson & Dreger?

[rmodrak]

For me the easiest way out without breaking the seiscomp interface (assuming it's tested on their side and actually gives correct results for the moment tenors) would be to add a new definition option which only differs in the sign of m1 and m4 (https://github.com/krischer/instaseis/blob/master/instaseis/database_interfaces/base_instaseis_db.py#L133-L139).

Certainly, this sounds like a good way forward.

But I would be more comfortable if I understood why this is the solution. Do you have another source for the definitions of the sources assumed by Minson & Dreger?

A better understanding here would also make me more comfortable. Regarding fundamental sources, Minson and Dregor 2008 cite Langston 1981, but that study doesn't actually provide details. I'll check with Carl to see if he could follow up with Minson.

The other reference I've looked at is the README and code from the FK package from Lupei Zhu. Because the implementation there is hard to follow, it's probably not the most useful reference for us, unfortunately.

[rmodrak]

I compared the Green’s functions used by instaseis and FK. All the Green’s functions agreed, up to the sign differences given in the table below. The instaseis and FK Green’s function definitions are given relative to the definitions required to make the MinsonDreger2008 formulas work, assuming MinsonDreger2008 uses a north,east,down basis, which they are not entirely clear about. Understanding the FK source code was difficult, and it is still necessary to double check the table.

The bottom line though is that every reference we’ve checked seems to use a different definition for the fundamental sources. This lack of standardization (or indeed absence of any definition from many studies, including MinsonDreger2008 or Langston1981) seems to be prevalent.

    instaseis       FK   MinsonDreger2008
ZSS     1           -1           1
ZDS    -1           -1           1
ZDD     1            1           1
ZEP     1            1           1

RSS     1           -1           1
RDS    -1           -1           1
RDD     1            1           1
REP     1            1           1

TSS    -1           -1           1
TDS     1           -1           1

[rmodrak]

Email exchange with Carl [edited for length]:

Hi Ryan,

I agree that this is a mess, both in the literature and in people's codes.

CAP attempts to explain the formulas that it uses, but I also don't see a clear match with the published equations.

In any case, this may not be resolvable by us. Are you comfortable (enough) using whatever tweak is needed for instaseis to match capuaf, then proceeding for now?

If you think we should engage Lupei Zhu, let me know. He would be the one to ask about cap and fk. I noticed a new function is his 2015 CAP release that may be relevant.

Thanks, Carl

Hi Carl,

Your phrase "using whatever tweak is needed to match capuaf" exactly describes the current state of the mtuq code. I was considering going just slightly further, and submitting a new get_greens_function pull request for instaseis. This would result in an mtuq speedup of about 10 percent, because we'd no longer need to perform ad hoc permutations and sign changes inside the grid search loop. How does this sound?

If you like, I am comfortable going ahead with the above plan to submit a relatively short pull request to instaseis, and leaving it there. I feel we have pursued this just about as far as the current literature allows, and further correspondence with Lupei Zhu or others may have diminishing returns.

Thanks, Ryan

[martinvandriel]

I am comfortable going ahead with the above plan to submit a relatively short pull request to instaseis, and leaving it there.

Same here, I think we digged deep enough to move forward.