ThomasLecocq/msnoise-tomo

Looking for guidance on SEEDing with the C++ code

Opened this issue · 2 comments

I have been working to modify iftan.py and ftan.py to work with any of the combination (PV, FV, PT, FT). I have everything working now, but the velocity seeding is broken. I am comparing my updated code to the version in the WHL file. Below is the iftan output for a SAC file.

fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.5
disp: no
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
write amplitude
write phase
write FT
write TV
fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.5
disp: cont,  finit=13.2006 vginit=0.185903
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
freq near 13.2006 = 13.2007
write amplitude
write phase
write FT
write TV

This is the output from running ftan() twice. ftan() is the external C routine. The first time the ftan() function is ran, you can see that the the "disp" is set to "no". The second time, the ftan() function is ran with the seed "finit" and vginit". The ftan() function is actually invoke in ftan_call.py. See below

ftan(filename, fmin, fmax, vgmin, vgmax, bmin, bmax,
         diagramtype, nfreq, ampmin, dist, disp="cont", tinit=pinit, vginit=vinit)

Now, in the modified version of my code, which is compiled on my mac, I cannot get the seeding to work properly. Here is the output for the exact same SAC file and db parameters. You can see I have added a couple of outputs, and a print of the vinit value before the second call to ftan().

Running FTAN the first time: /Users/dmikesell/Desktop/test_for_john/msnoise_geopark/TOMO_SAC/01/ZZ/GP_01101_GP_01205_MEAN.sac

fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.35
disp: no
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
write amplitude
write phase
write FT
write TV
Seed velocity: 0.185903 [km/s]
Seed period: 0.075754 [s]
vinit: 0.185903

Running FTAN a second time with the SEED

fmin=7, fmax=30, nfreqs=40
bmin=0.02, bmax=2
ampMin=0.3
vgMin=0.07 vgMax=0.35
disp: cont,  finit=13.2006 tginit=2.24583e-314
diag: Period - Vgroup
out: enable (matrix)
max: 13.2007
freq near 13.2006 = 13.2007
write amplitude
write phase
write FT
write TV

I am using the same diagramtype (i.e. PV), but the C code always spits out tginit=0, rather than vinit, which is what I am passing to it. My call is identical (I think) to that in the WHL file.

ftan(filename, fmin, fmax, vgmin, vgmax, bmin, bmax,
     diagramtype, nfreq, ampmin, dist, disp="cont", tinit=pinit, vginit=vinit)

However, the output is not the same. How can I tell which version of the code was used to create the WHL file? Or can someone tell me why the C code keeps looking for tginit (the time seed) instead of the velocity seed (vginit) when I use the PV diagram type?

I should add, I am running the whl file in py37 and my modified code in py38, both with numpy 1.18.1. I will do some tests of my modified code in the py37 environment and see if that has an effect.

This is not a python 3.7 vs. 3.8 issue. I get the same behavior on both. I cannot tell if it is my c compiler or not. I do not get any strange errors during compiling with gcc 8.3.0.