180 degree difference between SAzim and RAz

[TomWinder]

My understanding of the documentation is that SAzim and RAz should be identical when using travel-time grids constructed from 1D velocity models. Instead (in both NLLoc 6.0.0 and 7.0.0 - I haven't cloned this latest version) I find that RAz is 180 degrees different (i.e. opposite). This is the case irrespective of RDip.

Is this somehow intentional? It seems to me that the station azimuth calculation is correct, and the ray azimuth is wrong by 180 degrees.

I noticed this when parsing hyp files into MTfit - it turns out one of my predecessors in our research group had written a patch for NLLoc v6 to "correct" this, but when I switched to the vanilla version of NLLoc v7 the problem re-appeared (and all my moment tensors were upside down). I have attached the patched version of GridLib.c here, but to my (admittedly untrained) eye it looks quite clunky, so I'm not sure it will be of much use.

The Vel2Grid / Grid2Time parameters I have used are as follows:

TRANS LAMBERT WGS-84 63.8 -19.0 64.50 65.00 0.0

VGGRID 2 1000 400 0.0 0.0 -4.0 0.5 0.5 0.5 SLOW_LEN

LAYER -3.0 1.500 0.625 0.843 0.3511 0.0 0 LAYER 1.0 4.000 0.600 2.247 0.3371 0.0 0 LAYER 3.0 5.200 0.450 2.921 0.2528 0.0 0 LAYER 5.0 6.100 0.304 3.427 0.1708 0.0 0 LAYER 6.0 6.404 0.179 3.598 0.1006 0.0 0 LAYER 7.0 6.583 0.105 3.698 0.0590 0.0 0 LAYER 8.0 6.688 0.104 3.757 0.0584 0.0 0 LAYER 9.0 6.792 0.104 3.816 0.0584 0.0 0 LAYER 10.0 6.896 0.047 3.874 0.0266 0.0 0 LAYER 13.0 7.038 0.019 3.954 0.0107 0.0 0 LAYER 14.0 7.057 0.019 3.965 0.0107 0.0 0 LAYER 15.0 7.076 0.019 3.975 0.0107 0.0 0 LAYER 16.0 7.095 0.019 3.986 0.0107 0.0 0 LAYER 17.0 7.114 0.019 3.997 0.0107 0.0 0 LAYER 18.0 7.133 0.019 4.007 0.0108 0.0 0 LAYER 22.0 7.210 0.019 4.051 0.0107 0.0 0 LAYER 26.0 7.286 0.019 4.093 0.0107 0.0 0 LAYER 30.0 7.362 0.009 4.136 0.0053 0.0 0 LAYER 42.0 7.400 0.000 4.157 0.0000 0.0 0

GTMODE GRID2D ANGLES_YES

GTSRCE HOLR LATLON 64.87495 -16.81357 -0.713 0 (etc.)

GT_PLFD 1.0e-3 0

[alomax]

Hi Tom

The SAzim is measured for a line from the epicenter to the station, and the RAz is measured for a ray leaving the hypocenter towards the station. So, as you note, they should be the same for a 1D, laterally homogeneous model,

SAzim (float*6.2)
    Maximum likelihood or expectation hypocenter to station epicentral azimuth in degrees CW from North 
RAz (float*5.1)
    Ray take-off azimuth at maximum likelihood hypocenter or expectation in degrees CW from North 

I checked some NLL (NLLoc:v7.00.10(30Aug2020)) output I am running today for a 1D model and these two values are the same (though printed to different decimal precision).

Can you check or send examples of your event *.hyp output to see if this is the case? (I suppose you already checked this, but I wanted to confirm...)

Best regards,

Anthony

[TomWinder]

Hi Anthony,

Thanks for your quick response - I'm glad that I understood correctly. Yes of course - please find an example attached. I have found this consistently for a few 1D velocity models in different regions, though admittedly the scripts for generating the LUT's and running NLLoc itself are all closely related.

I have also attached the angle and time LUT files for one station to see if you can re-create the issue from those.

Tom

loc.20140829.201636.grid0.loc.hyp.zip HOLR_LUT.zip

[alomax]

Did you run Grid2Time and NLLoc on the same machine, with the same byte-order? I suppose so, otherwise there might be all sorts of other problems. And was Grid2Time run on the same version of NLL as NLLoc?

I cannot see a reason for the 180 deg flip, other than your rays are leaving the hypocenter away from the station, which seems unlikely. A way to check this is to plot the angles grid in Grid2GMT. You can check this, or I would need your control file to run Grid2GMT.

[TomWinder]

Yes, I re-made the grids with both v6 and v7, and ran then ran each with both NLLoc v6 and v7 to check it wasn't some sort of inter-compatibility issue. What I sent you here previously were grids created by Grid2Time v7 and the hyp file from running NLLoc v7 using those grids. All were made on the same machine.

I have attached my control file here (the reason I didn't before is due to the clutter.. sorry about that! but I have left it as-used in case there's something there that's causing the problem) and will have a go at plotting up the angles grids as soon as I get a moment.

dyke_tip.in.zip

[alomax]

I guess I need all the files (obs etc.) so I can run from scratch here (tomorrow...) and see if I get the same, and then debug... You can send by email alomax@free.fr

[TomWinder]

Will do. Thanks again for your help.

[alomax]

Well, when I run from scratch here (Vel2Grid, Grid2Time, NLLoc) I get the same angles for SAzim and RAz!

PHASE ID Ins Cmp On Pha FM Date HrMn Sec Err ErrMag Coda Amp Per > TTpred Res Weight StaLoc(X Y Z) SDist SAzim RAz RDip RQual Tcorr HOLR ? ? ? P U 20140829 2016 36.1800 GAU 2.00e-02 -1.00e+00 -1.00e+00 -1.00e+00 > 1.4403 -0.0431 1.5184 103.6066 121.6300 -0.7130 1.4028 139.05 139.1 160.8 9 0.0000 HOLR ? ? ? S ? 20140829 2016 37.2100 GAU 2.00e-02 -1.00e+00 -1.00e+00 -1.00e+00 > 2.5637 -0.1365 1.1496 103.6066 121.6300 -0.7130 1.4028 139.05 139.1 160.8 9 0.0000 RIMA ? ? ? P U 20140829 2016 36.3500 GAU 2.00e-02 -1.00e+00 -1.00e+00 -1.00e+00 > 1.5791 -0.0119 1.6929 100.1212 124.5224 -0.7460 3.1534 305.54 305.5 145.2 9 0.0000 RIMA ? ? ? S ? 20140829 2016 37.4400 GAU 2.00e-02 -1.00e+00 -1.00e+00 -1.00e+00 > 2.8109 -0.1537 1.0962 100.1212 124.5224 -0.7460 3.1534 305.54 305.5 145.2 9 0.0000 ...

What CPU and OS are you using? Which compiler?

Perhaps try the version here: https://github.com/alomax/NonLinLoc

Thanks, Anthony

[TomWinder]

How intriguing... I am compiling / running on Ubuntu 16.04 with gcc 5.4.0:

gcc --version
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.12) 5.4.0 20160609
Copyright (C) 2015 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

CPU is:

cat /proc/cpuinfo | grep "model name"
model name      : Intel(R) Core(TM) i5-7600 CPU @ 3.50GHz

I have cloned this repository. Which Makefile should I use to compile it? When I downloaded the source code for 7.0.0 from your website there was just a single Makefile, but here there are 3 - it looks like I should use either Makefile-NLLoc or Makefile.ORIG to compile with gcc.

[alomax]

There are instructions using cmake on the home github page.

[TomWinder]

I updated to cmake 3.18.2:

cmake --version
cmake version 3.18.2

And followed the build instructions, but unfortunately make failed with an error. Log files attached here.

[alomax]

Looks like it only failed for Loc2ddct, which is not needed for the main NLL tools. Are the other executables fully built? I know little about cmake, but it is supposed to be smart and help with building on diverse platforms...

[TomWinder]

Yes, Vel2Grid, Grid2Time and NLLoc were all there, so I have run through again from scratch. It still gives me the 180 degree difference: v7.0.1_loc.20140829.201636.grid0.loc.hyp.zip

For completeness, the list of binaries available was as follows (though I imagine it is not relevant here):

tebw2@cucurbita:/space/tebw2/Software/NonLinLoc/src$ ls bin/
fmm2grid  grid2scat  Grid2Time  GridCascadingDecimate  NLLoc  PhsAssoc  scat2latlon  sphfd_SWR_NLL  Time2Angles  Vel2Grid

[alomax]

OK. I am stumped. Could it be a subtle machine architecture or compiler issue? I think the only place where 180 deg is subtracted from the azimuth is in GridLib.c: /* determine azimuth (2D grids) */ if (gdesc.type == GRID_ANGLE_2D) { if (*pazim > 0.0) *pazim = sta_azim; else { *pazim = sta_azim - 180.0; if (*pazim < 0.0) *pazim += 360.0; } } You could try reversing these tests in the C source. These tests of *pazim are needed as to allow the case of the ray leaving the hypocenter away from the station for 2D time/angle grids.

[TomWinder]

I made the following change:

    /* determine azimuth (2D grids) */
    if (gdesc.type == GRID_ANGLE_2D) {
        if (*pazim < 0.0)
            *pazim = sta_azim;
        else {
            *pazim = sta_azim - 180.0;
            if (*pazim < 0.0)
                *pazim += 360.0;
        }
    }

i.e. reversing the test for *pazim > 0.0 (hopefully this is what you meant!), re-compiled and re-ran both from the previous grids and from scratch, all with no effect on the end result. I also can't immediately think why the ray would ever leave the hypocentre away from the station for a 1D model, so don't quite understand the test (but I am not particularly familiar with how the grids are set up and handled).

I'm afraid I don't have any bright ideas, other than perhaps compiling using a more up-to-date version of gcc, which I will try now. Did you use exactly the same control file when you ran it on your machine? And did you compare the angles grids I sent to those generated when you re-ran it?

I had a brief look into the Grid2GMT docs - am I correct in thinking that they require a 3D time grid to work? Either way I am not 100% sure I understand how I would plot the angles grid.

[TomWinder]

Aha - good news - updating gcc --> v7.3 seems to have fixed the issue! I reverted the source code before re-compiling. Full info on compilers used:

tebw2@cucurbita:/space/tebw2/Software/NonLinLoc/src$ cmake .
-- The C compiler identification is GNU 7.3.0
-- The CXX compiler identification is GNU 5.4.0
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working C compiler: /home/tebw2/.miniconda3/envs/qmigrate/bin/x86_64-conda_cos6-linux-gnu-cc - skipped
-- Detecting C compile features
-- Detecting C compile features - done
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++ - skipped
-- Detecting CXX compile features
-- Detecting CXX compile features - done
-- Configuring done
-- Generating done
-- Build files have been written to: /space/tebw2/Software/NonLinLoc/src

I also checked that running the newly compiled NLLoc using the old (vanilla v7.0.1 compiled with gcc 5.4) LUT still does have the 180 issue, so it seems to be in the generation of the angles grids.

I have attached the new, correct hyp file here.

[alomax]

Good, and thanks for all the info. Bizarre - maybe related to the way NLL compresses azimuth, dip and ray-angle quality into an integer that gives problems, otherwise I do not think there is anything that should cause hiccups at the compiler/machine level...

Just give an update that all is fine, after you have done more work...

Thanks,

Anthony

[TomWinder]

I did a little more digging, and things took a turn for the slightly more bizarre: I returned to the NLL v7.0.0 source-code and compiled that with gcc v7.3 and in this case it did not fix the issue.

So I un-commented the line: printf("DEBUG: ReadTakeOffAnglesFile: x y z %f %f %f raz %f rdip %f rq %d\n", xloc, yloc, zloc, *pazim, *pdip, *piqual); in GridLib.c to see if that would shed more light on what was going wrong. It shows that for the configuration where it "just works" (= NLL v7.0.1 compiled with cmake / gcc v7.3) the raz returned by GetTakeOffAngles(&angles, pazim, pdip, piqual) is positive (appears to always be the same, and nonsensical from my point of view):

DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 1.402803 5.620715  raz 6192.700000  rdip 160.800000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 4.202287 5.620715  raz 6192.700000  rdip 136.200000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 3.153422 5.620715  raz 6192.700000  rdip 145.200000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 3.757604 5.620715  raz 6192.700000  rdip 140.000000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 4.334428 5.620715  raz 6192.700000  rdip 135.500000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 4.353085 5.620715  raz 6192.700000  rdip 134.900000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 6.081384 5.620715  raz 6192.700000  rdip 119.700000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 6.396265 5.620715  raz 6192.700000  rdip 117.000000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 6.462429 5.620715  raz 6192.700000  rdip 117.100000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 7.083021 5.620715  raz 6192.700000  rdip 113.000000 rq  9

Whereas for all the configurations (combinations of compiler, version etc.) where it doesn't work, it returns raz == 0.000; e.g.:

DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 1.402803 5.620715  raz 0.000000  rdip 160.800000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 4.202287 5.620715  raz 0.000000  rdip 136.200000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 3.153422 5.620715  raz 0.000000  rdip 145.200000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 3.757604 5.620715  raz 0.000000  rdip 140.000000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 4.334428 5.620715  raz 0.000000  rdip 135.500000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 4.353085 5.620715  raz 0.000000  rdip 134.900000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 6.081384 5.620715  raz 0.000000  rdip 119.700000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 6.396265 5.620715  raz 0.000000  rdip 117.000000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 6.462429 5.620715  raz 0.000000  rdip 117.100000 rq  9
DEBUG: ReadTakeOffAnglesFile: x y z 0.000000 7.083021 5.620715  raz 0.000000  rdip 113.000000 rq  9

This of course explains why me reversing the test to:

    /* determine azimuth (2D grids) */
    if (gdesc.type == GRID_ANGLE_2D) {
        if (*pazim < 0.0)
            *pazim = sta_azim;
        else {
            *pazim = sta_azim - 180.0;
            if (*pazim < 0.0)
                *pazim += 360.0;
        }
    }

didn't make a difference, as the case that *pazim == 0.0 slips through with both > and < as the test.

So the question is, should this inequality really be "greater than or equal to", as opposed to "greater than"? For the very short term I have switched this to "greater than or equal to" and all versions compiled with all compiler versions now seem to work (including inter-compatibility between NLL v6.0.0 LUT's and NLL 7.0.0 / 7.0.1 NLLoc). But I would like to hear your thoughts on:

1. whether this is a valid solution
2. whether the values being returned by GetTakeOffAngles() - either 0 or a very large value, unrelated to the actual azimuth as far as I can tell - indicate that on my machine (with either compiler version) there is an error in the routines for compression / reading from an integer value; I would be interested to know what value of raz is returned when you run through on your machine (and/or to hear some more context for what this value actually means).

As a final note, I had a look at the test script included in this repository to see if it caught this bug: it does, but as currently written the script does not run. I had to make some fixes (adding lines to make directories for the time and model files, etc.), and I think the files it said to diff between were not quite correct. If it would be useful, I can create a separate issue to explain these problems in full, or just submit a PR with the fixes - whichever you would prefer.

Overall this episode has highlighted to me the importance of tests to catch this kind of subtle (but important!) bug!!

[alomax]

Hi Tom

I took a good look today, and the basic problem is deeper, there are 2 bugs that cancel each other, but with your compiler, one of the bugs changes an (invalid) unsigned short value and thus the results of the "if (*pazim" test we were looking at above. So I have made changes in several places to fix both bugs, but would like to do more tests here to confirm. The problem would affect the case where the ray leaves the hypocenter away from the station with any compiler, and also in the normal case where the ray leaves the hypocenter towards the station with your problem compiler/version. The bugs only affect 2D time grids (1D velocity models), not 3D models/time grids. And the 1D model would need to be very weird velocity model, so the rays leave the hypocenter away from the station (not really possible with a 1D, flat layered model).

I think you can most simply make a work around by changing: in GridLib.c around line 6336

        //dip = atan2(grady, -gradz) / cRPD;
        dip = atan2(fabs(grady), -gradz) / cRPD;

as long as you do not have a weird velocity model .

Tell me if this works (provisionally).

Anthony