diff --git a/src/config.yaml b/src/config.yaml index 3a8974e..48e676a 100644 --- a/src/config.yaml +++ b/src/config.yaml @@ -297,13 +297,16 @@ Inversion: # to be performed. Multiple rules can be defined simultaneously. # The best-correlated solution will be chosen, however, # more time will be needed for finalizing the MT computation. - # For example, if 4.0<= Mag <=4.5, Gisola will use 245.76 sec as + # For example, if 4.0<= Mag <=4.5, Gisola will use 204.8 sec as # inversion time window. Values, that can be used must be able - # to be divided with 8192 and lead to not repeating decimals. + # to be divided with 1024 and lead to not repeating decimals. Window: - - [4.0, 4.5, 245.76] - - [4.6, 5.5, 327.68] - - [5.6, 9.0, 409.6] + - [3.0, 3.4, 102.4] + - [3.5, 3.8, 102.4] + - [3.8, 4.5, 307.2] + - [4.6, 5.5, 409.6] + - [5.6, 7.1, 409.6] + - [7.2, 9.0, 1024] # Here is defined the time search dimension in the MT search grid. # Multiple rules can be defined simultaneously. The best-correlated @@ -315,16 +318,16 @@ Inversion: # relative to initial origin's time. # To determine the exact time in sec you have to associate it with the # corresponding values from Window parameter. For example, - # if 4.0<= Mag <=4.5, Gisola will use 245.76 sec as inversion time window, - # leading to a time kvantum of 246.76/8192=0.03 sec. Thus, -67*0.03=-2.01 sec - # to 5.01 sec with a step of 0.3 sec. A total of 24 time points for each 3D grid point. + # if 4.0<= Mag <=4.5, Gisola will use 204.8 sec as inversion time window, + # leading to a time kvantum of 204.8/1024=0.2 sec. Thus, -67*0.2=-13.4 sec + # to 33.4 sec with a step of 2 sec. A total of 24 time points for each 3D grid point. # In our case would lead to 24*1271 grid points= 30504 MT inversions. # Time Shifts must also comply with the parameters defined in the # parameters.f90 file within the inversion code. TimeShift: - - [4.0, 4.5, [-67, 10, 167]] # must be the same as the one in parameters.f90 - - [4.6, 5.5, [-51, 10, 161]] - - [5.6, 9.0, [-81, 15, 341]] + - [3.0, 4.5, [-10, 1, 20]] # must be the same with the one in parameters.f90 + - [4.6, 5.5, [-5, 1, 16]] + - [5.6, 9.0, [-5, 2, 43]] # Here is defined the frequency band that will be used for the inversion. # In form of F1 F2 F3 F4, where F1-F2 and F3-F4 are tapered. diff --git a/src/core/inversion/isola15.f90 b/src/core/inversion/isola15.f90 index f6d9fcb..2522044 100644 --- a/src/core/inversion/isola15.f90 +++ b/src/core/inversion/isola15.f90 @@ -209,7 +209,7 @@ PROGRAM ISOLA15 93 close(110) KEYINV=2 !fixed to deviatoric type - DT=TL/8192. + DT=TL/REAL(NUM_OF_TIME_SAMPLES) ISOURMAX=NS ISUBMAX=1 VARDAT=2.0e-12 diff --git a/src/core/inversion/isola15.f90.old b/src/core/inversion/isola15.f90.old index c2dc975..ccd9066 100644 --- a/src/core/inversion/isola15.f90.old +++ b/src/core/inversion/isola15.f90.old @@ -209,7 +209,7 @@ 93 close(110) KEYINV=2 !fixed to deviatoric type - DT=TL/8192. + DT=TL/REAL(NUM_OF_TIME_SAMPLES) ISOURMAX=NS ISUBMAX=1 VARDAT=2.0e-12 diff --git a/src/core/inversion/parameters.f90 b/src/core/inversion/parameters.f90 index 3ef3685..fac9d31 100644 --- a/src/core/inversion/parameters.f90 +++ b/src/core/inversion/parameters.f90 @@ -4,7 +4,7 @@ MODULE PARAMETERS INTEGER, PARAMETER :: MAX_SOURCE_POSITIONS = 100 - INTEGER, PARAMETER :: NUM_OF_TIME_SAMPLES_LOG2 = 13 + INTEGER, PARAMETER :: NUM_OF_TIME_SAMPLES_LOG2 = 10 INTEGER, PARAMETER :: NUM_OF_TIME_SAMPLES = 2**NUM_OF_TIME_SAMPLES_LOG2 diff --git a/src/isola.py b/src/isola.py index 5bca1de..a6ff802 100644 --- a/src/isola.py +++ b/src/isola.py @@ -37,6 +37,8 @@ # local import config, event, modules.kagan +NUM_OF_TIME_SAMPLES = 1024 + def kaganCalc(x): return modules.kagan.get_kagan_angle(*x[0:3], *x[3:6]) @@ -102,7 +104,7 @@ def calculateVariance(observed, synthetic, tl): """ with np.errstate(divide='raise'): try: - dt = round((tl/8192.0),4) + dt = round((tl/float(NUM_OF_TIME_SAMPLES)),4) obs = np.array(observed) syn = np.array(synthetic) @@ -340,11 +342,11 @@ def createRaw(): os.makedirs(os.path.join(config.rawdir,str(i))) _st=config.st.copy() - # downsampling at known frequency at 8192 elements - _st.resample(8192/tl) + # downsampling at known frequency at NUM_OF_TIME_SAMPLES elements + _st.resample(NUM_OF_TIME_SAMPLES/tl) # be sure than no point exceeds number for tr in _st: - tr.data=tr.data[:8192] + tr.data=tr.data[:NUM_OF_TIME_SAMPLES] for stat in stations: _st2=_st.select(station=stat.split()[4]) @@ -359,7 +361,7 @@ def createRaw(): ) else: # dummy data - for _ in range(8192): + for _ in range(NUM_OF_TIME_SAMPLES): text+='{:.6e}\t{:.6e}\t{:.6e}\t{:.6e}\n'.format(0, 0, 0, 0) with open(os.path.join(config.rawdir,str(i),stat.split()[4].upper()+'raw.dat'), 'w') as _: @@ -564,7 +566,7 @@ def gatherResults(cfg=None, evt=None, workdir=None, bestinvdir=None, revise=Fals org=Origin() mag=Magnitude() - org.time=event.getOrigin(config.cfg,evt,config.cfg['Watcher']['Historical']).time+float(best[4])*(tl/8192) + org.time=event.getOrigin(config.cfg,evt,config.cfg['Watcher']['Historical']).time+float(best[4])*(tl/NUM_OF_TIME_SAMPLES) _lat, _lon=cart2earth(event.getOrigin(config.cfg,evt,config.cfg['Watcher']['Historical']).latitude, event.getOrigin(config.cfg,evt,config.cfg['Watcher']['Historical']).longitude, best[1], best[2]) org.latitude=_lat diff --git a/src/plot.py b/src/plot.py index 97d8eba..b357329 100644 --- a/src/plot.py +++ b/src/plot.py @@ -44,6 +44,8 @@ from obspy.core.event import read_events import multiprocessing +NUM_OF_TIME_SAMPLES = 1024 + def emsc(filepath='emsc.txt'): filepath=os.path.join(config.outputdir, ('emsc.txt' if not config.revise else 'emsc.revise.txt')) @@ -564,7 +566,7 @@ def top(solutionspath=None, tl=None, filepath='top.png'): '.hed'),'r') as f: tl=float(f.readlines()[2].split('=')[1]) - dt=tl/8192 + dt=tl/NUM_OF_TIME_SAMPLES if solutionspath: with open(solutionspath, 'r') as f: solutions=f.readlines() @@ -692,7 +694,7 @@ def major_formatter(x, pos): '.hed'),'r') as f: tl=float(f.readlines()[2].split('=')[1]) - dt=tl/8192 + dt=tl/NUM_OF_TIME_SAMPLES if solutionspath: with open(solutionspath, 'r') as f: diff --git a/test/benchmark/inventory.xml b/test/benchmark/inventory.xml index 1364bae..4772c95 100644 --- a/test/benchmark/inventory.xml +++ b/test/benchmark/inventory.xml @@ -17361,7 +17361,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -17678,7 +17678,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -17995,7 +17995,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -24882,7 +24882,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -25199,7 +25199,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -25516,7 +25516,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -25793,7 +25793,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -26070,7 +26070,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -26347,7 +26347,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -36220,7 +36220,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -36538,7 +36538,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -36856,7 +36856,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -38945,7 +38945,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -39222,7 +39222,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -39499,7 +39499,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -39776,7 +39776,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -40053,7 +40053,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -40330,7 +40330,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -48109,7 +48109,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -48427,7 +48427,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -48745,7 +48745,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -49022,7 +49022,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -49299,7 +49299,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -49576,7 +49576,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -51689,7 +51689,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -51987,7 +51987,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 @@ -52285,7 +52285,7 @@ COUNTS ODD - 8192.0 + 1024.0 98304.0 540672.0 1802240.0 diff --git a/test/config.yaml b/test/config.yaml index 38304ae..bfc823d 100644 --- a/test/config.yaml +++ b/test/config.yaml @@ -28,7 +28,7 @@ Inventory: Distance: - [6.1, 9.0, [100, 700, ['BH','HH']]] - Azimuth: [3,2] # sectors, stations per sector # null for all! + Azimuth: [2,2] # sectors, stations per sector # null for all! # The Stream Service supports FDSNWS-station only and file in XML format (as mentioned in ObsPy documentation) Stream: @@ -62,9 +62,12 @@ Green: Inversion: - + # InversionWindow = 409.6, + # dt = 409.6/1024 = 0.4 sec + # TimeShift: -10*0.4 = -4 ---- 43*0.4 = +17.2sec + # Step: 0.4 * 2 = 0.8 TimeShift: - - [5.6, 9.0, [-81, 15, 341]] + - [5.6, 9.0, [-10, 2, 43]] Window: - [5.6, 9.0, 409.6]