Throughout my PhD work, I’ve built several pieces of research code, mostly in c. I feel writing software is an integral part of any science PhD, because if the necessary software already exists, then the problem is solved and you should work on something else. But also it really helps you think through the problem and better understand why things work or don’t work.
Matlab R2014b has implemented a new graphics program. This has led to a few bits to not work quite properly in SplitLab and FuncLab. Download these versions to avoid these problems: [SplitLab1.2.1.zip] [Funclab1.8.2.zip] Note that these have not been fully tested, so bugs may still exist. FuncLab is also available with topography and coastline data for use with the Mapping Toolbox [Funclab1.8.2_with_map_data.zip].
Along these lines, Matlab R2014b comes with a new default color palette called “Parula”. This is a really fantastic color palette, and if you want to use it in GMT you can download it here: parula.cpt
FuncLab I am working on updating FuncLab to improve functionality and added customization. Interface updates include color palette swapping and new plot tools to look at stacked Receiver Functions. Work in the background include improved version tracking within a given project, depth mapping based on a 1D velocity model, station and event maps without the Matlab Mapping Toolbox, Receiver Function computations within FuncLab, and much more. For a video tutorial of the new FuncLab, check out this YouTube video from the IRIS EPO. The manual has recently been updated for version 1.8.0 and is available here [funclab_manual_1.8.0.pdf]. This project is not yet complete, so some of the new functions are not yet available. [FuncLab1.8.2.zip]. FuncLab is also available with topography and coastline data for use with the Mapping Toolbox [Funclab1.8.2_with_map_data.zip].
If you use FuncLab, please cite our paper: Porritt, R. W. and Miller, M. S., (2018), Updates to FuncLab, a Matlab based GUI for handling receiver functions. Computers and Geoscience, 111, 260-271, doi:10.1016/j.cageo.2017.11.022
Michigan receiver function data example: [here]
Station Analysis Tools One useful pack I wrote is a set of c routines for computation of power spectral densities, coherence, probability density functions, and a handful of other tools for monitoring the health of a station. These are effectively like PQLX or IRIS’s Quack, but are designed for very simple inline processing with sac data files and associated pole-zero response files. I’ve packaged these tools with the standard unix method, so just download the tarball, untar, configure, and make and you should be ready to go. Includes example data files, but requires the proper installation of fftw3 and pkgconfig. [station_analysis_tools]
[Old] SplitLab1.2.1 SplitLab is the standard shear-wave splitting analysis tool built in Matlab. It was written in 2006 by Andreas Wüstefeld and was due for an update. The original program and a splitting analysis database can be found here: http://splitting.gm.univ-montp2.fr/. I’ve updated it to utilize irisFetch.m for obtaining event and waveform data, removed the dependence on the mapping toolbox, and created a new output called ‘station.mat’ which contains the full misfit space for each splitting measurement. Download here: [SplitLab1.2.1.zip]. Please note that I’ve stopped working on this code as the original French group have been making strong strides in improving SplitLab.
Rob Porritt 
Hi Rob. Would you like us to put bug reports here?
Hi Shahar. This will work for now. Make sure to include the red error information if possible. This medium allows myself and other users to discuss the codes.
SplitLab: I suggest to add the following line to SL_changeServiceDataCenter:
case ‘http://www.data.scec.org’
config.serviceDataCenter = ‘http://service.scedc.caltech.edu’;
And perhaps even change “www.data.scec.org” to “scedc.caltech.edu” in the calling function, since data.scec.org does not exist anymore.
Thanks! Updates made and version 1.2.1 uploaded.
Your current Matlab version 8.5.0.197613 (R2015a)
———————————————————————-
Funclab 1.7.7
Adding FuncLab1.7.7 receiver function seismic tools to your MATLAB paths
You have the Mapping, Signal Processing Toolbox, and Parallel computing toolboxes, which provide full
access to all the functionality of the FuncLab package.
Error using textscan
Invalid file identifier. Use fopen to generate a valid file
identifier.
Error in setup_funclab (line 181)
javaClassPathCell = textscan(javaClasspathFid,’%s’);
Sorry about that; new version being uploaded that should fix this.
Mr.Shahar,did joint inversion is also the open source code ?
For joint inversions, I’d start with looking into Bob Herrmann’s Computer Programs for Seismology. He has good examples setup for you to run through and the codes are fairly easy to use.
Known bug with SplitLab1.2.1 and Matlab versions 2014b and newer: after associating events, if you have the mapping toolbox, the statistics plotting tool will fail due to a change in plotm and its interaction with matlab plotting primitives. To work around, force the lines which test if the mapping toolbox is available to evaluate to false.
I’m sorry, how do I “force the lines which test if the mapping toolbox is available to evaluate to false.” I’m not any better with MatLab than I used to be.
Thank you, quite useful.
I was poking around the SplitLab code to try and make it work, but my extensive knowledge of disp(‘Hello world’) wasn’t enough.
I am a student and have a very basic experience with Linux, C and coding. But i wanna use either PQLX or your station_analysis_tools software for PSD and PDF. Could you help me to make a step by step installation guide so i may be able to use your station_analysis_tools software package in CentOS 6.3 32 bit version.
Warm regards,
Abid Ali
Biggest thing is to make sure you have FFTW installed. Its free and easy and I highly recommend it. To use the tools, just try them out on the command line by going to the directory and trying a command like:
> ./compute_psd
and it’ll give you the usage. Basically you’re going to need SAC formatted files and the SAC Polezero response file.
I am quite sorry to say that i tried using command line by ‘compute_psd’ but i did not fully comprehend how to give the input .SAC file see the command below
>compute_psd [ KS.ADO2..BGE.D.2014.001.000000.SAC]
but there was no output file i got . If you have some ‘getting started ‘ manual or you can show how to do that job it would be highly appreciable.
Thanks for your cooperation !
Check out the README file or the seiscode wiki page. But simply, the way it works is you need to create a new text file with specific parameters. Then your first argument on the command line will be the file name of that text file.
Hi Rob. I have downloaded SplitLab 1.0.5 from website you offered and found that the ATD example project data has something difficult. The theoretical travel time caculated by taup toolbox doesnot match the actual SKS phase. I thought it may be my fault in the installing and setting berore the process. But I am not sure about it.Could you please do me a favour to solve it ?
Hi Ju,I’ve found that sometimes there is a problem with the depth definition – ie some events are in meters and taup assumes km for depth. This could cause a noticeable offset if the event is deep.
Rob, I checked the event data head file and confirmed that the event depth is written in km style. But still the problem exists. However, the difference between taup and SKS phase of ATD data is not so big but not regular. So I don’t know what causes the difference.
How many seconds off are the prediction and the phase? If small, it may just be heterogeneity along the raypath. It may also just be that the event has a weak SKS arrival. However, I can’t be sure as I haven’t yet checked the example data you’ve mentioned.
The time difference may be about 10 seconds. And sometimes the “A” and “F” given by the author is before the therotical SKS travel time and sometimes after it. I don’t know the reason. So I thought maybe it’s my wrong during the installing or setting before processing.
Sorry I don’t have more chance to look into this right now. You might want to double check the SAC manual to find out the meaning of those header variables and their relation to timing of the seismogram. Try adjusting them either in sac or matlab and then testing out Splitlab’s and Taup’s capabilities.
Hi Rob. I am a student and I am trying to work out funclab fetch traces and create project. I’ve tried several computers and different matlab version but that’s not the solution. Below is the bug in matlab:
Matlab: Unable to create directory C:\Program Files\MATLAB\R2009b\toolbox\lib\
Please add matTaup.jar and IRIS-WS jar file to your local path and update
your classpath.txt file.
funclab:??? Error using ==> taupTime at 117
Java exception occurred! Please check model name.
Error in ==> funclab>FL_new_fetch_project_setup_Callback at 3138
tt = taupTime(‘iasp91′,eventDepth,traceStartPhase,’deg’,gcarc);
??? Error while evaluating uicontrol Callback
Could you please do me a favour to solve it ?
Hi Chih-Wei Cheng, sorry about the lack of a response. I may have fixed this error in the newest release of funclab 1.7.8.
Hey Rob,
during searching for updates regarding SplitLab I found your extended version.
I think for all users it would be nice to add the fixed determination of the degrees of freedom for the SC error estimation published by Walsh et al. (2013).
Walsh, E., R. Arnold and M. K. Savage (2013), “Silver and Chan Revisited”, JGR
best
Michael
Additionally I found a potential error in function preSplit.m, lines 143-145, where the original T component is used to calculate the errors instead of the resulting T component after linearisation. Can you please cross-check this?
bests
Michael
Thanks for the comments Michael. I haven’t had much chance to look into this, but will try to look into soon. I’ve got some other ideas to improve SplitLab too, just been busy recently with other projects.
Dear Rab,
I am using Funlab 1.7.7 for H-k and CCP stackings beneath Thailand (Lat 8-20, Long 98-105). I am not sure that DNA-13.mat could be used in CCP stacking for an area beneath Thailand.
Regards,
Kasemsak Saetang
Dear Rob,
I am using Funclab 1.7.7 for H-k and CCP stacking in Thailand (Lat 8-20, Long 98-105). I am not sure that DNA-13 tomographic model could be used for Thailand.
Regards,
Kasemsak Saetang
Hi, I’d suggest either finding and adding a model more appropriate for Thailand, or use a global model, such as GyPSuM which is included.
I just completely created a model beneath Thailand by using CCPV2 with the GyPSuM tomographic model. I don’t know how can I find a method used in CCPV2 for writing a scientific paper. Please give me more information about CCPV2.
Sincerely,
Kasemsak Saetang
Hey Rob, I am using Funclab 1.7.10 for H-k stacking in Caspian Sea. The Matlab version that I use is 2016a. I get the error when I am trying to create the project:
Error using taupTime (line 117)
Java exception occurred! Please check model name.
Error in funclab>FL_new_fetch_project_setup_Callback (line 3244)
tt = taupTime(‘iasp91′,eventDepth,traceStartPhase,’deg’,gcarc);
Could you please send a version that works with 2016a or maybe an older version where this is not an issue?
Best,
Elshan
Hi Elshan,
The problem is with the install of MatTaup. As far as I’ve found, there is not an easy way to install java packages with an automated installer.
Hello Rob
Is this Funclab working under matlab for windows ? I study about receiver function in Java Island, Indonesia, right now. Could you send me the manual user for your Funclab version ? Thank you
Yusuf
Hi Yusuf, the code should work under Windows, however a few file system errors may occur. It was designed on mac and linux. There is no manual yet, but I’m working on making it a bit more self-documenting. For instance, version 1.8.0 (just uploaded) has a menu with seminal citations for reference.
Hello Rob,
We are a group of students who are now studying geophysics (especially seismology) at Academia Sinica Taiwan. We are also trying to learn/practice writing program for our work. Could you have a visit in our website?
Thank you and wish you a nice day.
Dear Rob,
Could you facilitate me a gmt script to plot probability density spectra functions (general_pdf output)? Many thanks
Say the output is tmp.prob and you’re using GMT5, then:
gmt makecpt -Cgray -I -D -Z > tmp.cpt
gmt psxy tmp.prob -R${xmin}/${xmax}/${ymin}/${ymax} -JX7i/7i -Bxa${xannot}f${xannot} -Bya${yannot}f${ytick} -BnSeW -Ss0.1i -Ctmp.cpt > tmp.prob.eps
Hi Rob,
I am trying to use CCP Stacking V2 using my own 3D model. However, when I run the Ray Tracing I get the following error:
Error using griddedInterpolant
Interpolation requires at least two sample points in each dimension.
Error in interp1 (line 161)
F = griddedInterpolant(X,V,method);
Error in ccpv2_do_ray_tracing_serial_PtoS (line 396)
interpEposition = interp1(deps{1},epos{1},InterpDepthVector,’linear’,’extrap’);
Error in ccpv2_raytrace>saveraytrace_params_and_go_Callback (line 287)
ccpv2_do_ray_tracing_serial_PtoS(CCPV2RayTraceParams);
The size of deps{1} and epos{1} vectors varies between ~1000 to 1 for different stations. As expected, I get the error above when their size is 1. I couldn’t figure out what deps{1} and epos{1} are. The code in “ccpv2_do_ray_tracing_serial_PtoS” is a bit cryptic. I was hoping you might know.
Thanks for the help!
Hello Rob,
Great program! I’ve found it very useful.
Out of curiosity, why does it save the data twice? Once in /RAWTRACES, and once in /RAWDATA/RAWTRACES. It seems like this is a waste of computing time and storage space.
Just curious!
The original version of FuncLab makes a copy of the data into RAWDATA/ from wherever you have the data already stored. That subroutine is called again at the end of the fetching, so that copy command is executed. More importantly, I’ve decided to keep this system to make sure the user has a clean copy of the data either as backup or for other usages. It shouldn’t cost much computation time, but the extra storage can be an issue.
I’m backing everything up over Google Drive, and it makes a considerable difference if I have 5,000 files as opposed to 2,500 files. Do you know if there’s an easy way to disable the double back up? If not I’m sure I can find a workaround. It would just be nice to be able to do it by altering the source code.
Funclab 1.8.0. – feature suggestion
==========================-=
Currently, the auto trace editing options can be applied either to a selected table or on all the tables. In the latter option, it seems to compare traces from different stations. It would be useful to have an option to do (in one click) all stations separately.
Rob, thanks again for the good work you are putting into Funclab!
Hey Rob,
Like your program.
I am running it Funclab1.8 on Windows machine and tried to make a new project “from exisiting RFs, but run into following error:
Error using dataread
Delimiter has bad \ constant.
Error in strread (line 48)
[varargout{1:nlhs}]=dataread(‘string’,varargin{:}); %#ok
Error in funclab>new_project_setup_Callback (line 1547)
DirName = strread(DataDirectory,’%s’,’delimiter’,filesep);
Error while evaluating UIControl Callback
I have preprocessed the files according to the Funclab Manual and used your file/folder name convention. Maybe you might know.
Thanks for your help!
Hello Rob,
In order to use funclab functions, I added the paths to the files in the MATLAB session. When I run the setup.m in Matlab, it showing:
Warning: Unable to save path to file ‘/usr/local/MATLAB/R2015b/toolbox/local/pathdef.m’. You can save
your path to a different location by calling SAVEPATH with an input argument that specifies the full
path. For MATLAB to use that path in future sessions, save the path to ‘pathdef.m’ in your MATLAB
startup folder.
> In savepath (line 169)
In setup_funclab (line 378)
Please help me to set up the Funclab as I am a beginner.
Hello Rob,
The Funclab program is very useful for the research of receiver function. But I have a small question. I am using the version of 1.8.1 and when I conduct the manual trace edit, there is something wrong if I put the button of select all or deselect all which is Attempt to reference field of non-structure arra. How can I fix these mistake. Thank you.
Thanks for finding that – its a bug that I need to find and fix.
Hi Yaohui, I was unable to replicate this on my setup. What version of Matlab and operating system are you using?
Hi, Rob, the matlab version I am using is R2011b and the system is Ubuntu 16.04. Will they be OK for the program? Thank you!
Ubuntu should be fine, but I think that Matlab version may be too old. Sorry about that.
OK, maybe that is the problem and I will change a new one to have a try. Thank you for your advice and your program. Best wishes to you!
Hi Rob. Excelent program! its very useful :). I love the visualization part =)
Just two questions:
1-. I have a lot of waveforms in SEISAN format. Which is the best way (in your opinion) to convert SEISAN -> SAC format so that FuncLab can read it correctly?
2-. When I try to use Manual Trace Edit with the test-data, I can not plot the RF. This error message appears:
Error while evaluating UIControl Callback
I tried to do somethings but nothing works. All the others tools work properly.
Thank you very much!
Hi Rob,
I’m currently using funclab1.8.1 to perform RF stacking. This is a nice software.
I got a small question regarding the data format of the 3D seismic tomography model used for the 3D iteration in the RF to depth mapping. I noticed that there are three 3D tomography models in the directory, TOMOMODELS. I did some research into the data and the code by myself. It seems that the velocity perturbations of DNA13 and GYPSUM are expressed in percent value, while the velocity perturbation of NWUS_P09b_1-2 is not in percent value, since the ModelS and ModelP values in NWUS_P09b_1-2.mat are about 100 times smaller than those in DNA13.mat and GYPSUM.mat. In the velocity iteration code, it seems that the calculation assumes that the values in 3D tomography are in percent (for DNA13 and GYPSUM). I’m a little confused on this. If the perturbation values in NWUS_P09b_1-2.mat are not in percent, I cannot find the place in the code to take care of this. I may understand things wrong, but can you give me a little description on the data format of the 3D tomography models? I may put into other models after I can figure out the data format.
Hi Shangxin,
The depth mapping code checks a variable in the TOMO model structure called .TomoType. If this is set to “Perturbation”, like DNA13 and GyPSuM, it converts from a percentage value to an absolute velocity assuming the selected 1D model.
Cheers,
Rob
Hi Rob,
Yes. I noticed the function of the char variable “TomoType”. The current three 3D tomography models (DNA13, GyPSuM, and NWUS_P09b_1-2) all set this variable to “Perturbation”. The focus of my question is about the data format of the NWUS_P09b_1-2. The ModelS and ModelP variables of DNA13 and GyPSuM are indeed in percentage value (consistent with the code), but the perturbation values in NWUS_P09b_1-2 seem not in percent. When I load it in Matlab, the values are far smaller (~100 times smaller) than those of the other two tomography models. This makes me think that NWUS_P09b_1-2 expresses the perturbation in the absolute deviation (e.g., 0.01) instead of percentage (e.g., 1%). That is my confusion considering that the source code assumes the perturbation in percentage. I may understand things wrong but can you give a little hint on this?
Cheers,
Shangxin
Your understanding is correct. The NWUS_P09b_1-2 model is in absolute deviation (km/s) whereas GyPSuM and DNA13 are in percent. Based on their methodology, I’m not sure how they could have gotten an absolute deviation, so in version 1.8.2 I’ve upped their model by 100x to be on scale with the percentage models.
The biggest effect that this will have is in the depth mapping and ccpv2 which are routines I built and assumed velocity in percent. In ccp and gccp, the original authors wrote and assume absolute deviation. I’ve adjusted ccp and gccp to account for the increased values in NWUS_P09b_1-2.mat in the current, 1.8.2 version.
This comment and fix went up late at night and I probably wasn’t thinking as clearly as I should have. I have a better fix in mind and will put that up as soon as possible. This will allow for models with absolute perturbation (km/s), percentage perturbation (%), or absolute velocity (km/s) to be handled separately in all relevant instances.
Hi Rob,
Thanks for clarifying the issue. I need to report this to others in our group so just for confirmation: From your first reply, now NWUS_P09b_1-2 has been increased by 100x, and the subroutines (including the new depth mapping and ccpv2 written by you, as well as the original ccp and gccp) in funclab 1.8.2 version all process the velocity perturbation in percent consistently, yes?
Your mention of the unit km/s of the absolute perturbation still makes me a little confused. In my mind, the absolute deviation should still be a relevant value like 0.01 (corresponding to 1% for percent value) with no unit. I’m not sure whether this is a typo or something else.
I also agree with your follow-up fix. It’s of course better to indicate these three data types in the char variable of the 3D tomography model and allow for the separate handling.
Best,
Shangxin
Hi Rob,
I wonder whether the current funclab 1.8.2 enables the output of the results of the single station RF stacking trace into a text file? I’d like to plot it with something else so I need the data to process separately. Or is their a function to pick up the peaks from the single station stacking trace, say the 410 and 660 conversion peaks and then output the depth of the peaks for each stations into a text file?
Are these functions available now or I should write some plugin into the funclab to achieve this?
Best,
Shangxin