Matlab/Octave tools for geophysical studiesFrançois Beauducel 
Model  Description  Parameters  Outputs  References 
 

Space  Source  Medium  
Mogi  Point source in elastic halfspace Approximation for sphere of radius a << f 
r  f, ΔV f, ΔP, a 
ν E, ν or μ 
displacements, tilt, strains at surface  [Anderson, 1936] [Mogi, 1958] 
= 1.5  
Sun  Pennyshaped crack in elastic halfspace Approximation for h/a >> 1 
r  h, a, ΔV h, a, ΔP 
E, ν 
displacements at surface  [Sun, 1969]  
Okada  Rectangular dislocation in elastic halfspace  x,y  d, δ, θ, W, L R, U_{SLIP}, U_{OPEN} 
ν  displacements, tilt, strains at surface  [Okada, 1985] 

= 1 (sill) = 0.75 (dyke) 

Okubo  Rectangular dislocation in elastic halfspace  x,y  d, δ, θ, W, L R, U_{SLIP}, U_{OPEN} 
ν, ρ, ρ'  gravity and elevation change at surface  [Okubo, 1992] 
r  Radial distance from source 

x  Xcoordinate from source 
y  Ycoordinate from source 
f, h  Focal depth of source 
a  Radius of source 
d  Depth of topedge fault 
W  Width of fault 
L  Length of fault 
δ  Dipangle of fault 
θ  Strikeangle of fault (azimut) 
ΔV  Volume variation of source 
ΔP  Pressure variation of source 
R  Rakeangle of dislocation 
U_{SLIP}  Slipdislocation (in Rakeangle direction) 
U_{OPEN}  Tensiledislocation (opening) 
μ, G  Rigidity of medium (Lamé's second parameter) 
λ  Lamé's first parameter 
E  Elasticity of medium (Young's modulus) 
ν  Poisson's ratio of medium 
ρ  Density of medium 
ρ'  Density of cavityfilling matter 
Some conversion formulas between elastic parameters λ, μ, E and ν:
 
 
 
 

The Mogi [1958] model calculates analytical solution for surface deformation due to a point source in an elastic halfspace. This model is widely used to simulate ground deformation produced by local perturbation like volcanic magma chamber. It computes displacements, tilt and strain in a polar space, due to a volume variation at depth [Anderson, 1936] or an isotropic pressure variation in a spherical source [Mogi, 1958]. The proposed Matlab script is a literal transcription of Mogi's simple equations, extended to nonisotropic medium (Poisson's ratio different from 0.25). All parameters can be vectorized. See help for syntax, and script comments for details. MOGI Mogi's model (point source in elastic halfspace). [Ur,Uz,τ,σr,σt] = MOGI(R,F,ΔV,ν) or MOGI(R,F,A,ΔP,E,ν) computes radial and vertical displacements Ur and Uz, ground tilt τ, radial and tangential strain σr and σt on surface, at a radial distance R from the top of the source due to a hydrostatic pressure inside a sphere of radius A at depth F, in a homogeneous, semiinfinite elastic body and approximation for A << (center of dilatation). Formula by Anderson [1936] and Mogi [1958]. MOGI(R,F,ΔV) and MOGI(R,F,A,μ,ΔP) are also allowed for compatibility (Mogi's original equation considers an isotropic material with Lamé constants equal, i.e., λ = μ, ν = 0.25). Input variables are: F: depth of the center of the sphere from the surface, ΔV: volumetric change of the sphere, A: radius of the sphere, ΔP: hydrostatic pressure change in the sphere, E: elasticity (Young's modulus), ν: Poisson's ratio, μ: rigidity (Lamé constant in case of isotropic material). Notes:  Equations are all vectorized, so variables R,F,V,A,μ and P can be vectors or ND matrix of the same size and any of them can be scalars; then outputs will be also vector or matrix.  Convention: Uz > 0 = UP, F is depth so in Z direction.  Units should be constistent, e.g.: R, F, A, Ur and Uz in m imply ΔV in m³; E, μ and P in Pa; τ in rad, σr, σt and ν dimensionless. Example for a 3D plot of exagerated deformed surface due to a 1bar overpressure in a 10cm radius sphere at 1m depth in rock: [x,y] = meshgrid(3:.1:3); [th,rho] = cart2pol(x,y); [ur,uz] = mogi(rho,1,0.1,1e5,10e9,0.25); [ux,uy] = pol2cart(th,ur); ps = 1e8; surf(x+ux*ps,y+uy*ps,uz*ps), axis equal, light Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 1997 Updated: 20120405 References: Anderson, E.M., Dynamics of the formation of conesheets, ringdikes, and cauldronsubsidences, Proc. R. Soc. Edinburgh, 56, 128157, 1936. Mogi, K., Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them, Bull. Earthquake Res. Inst. Univ. Tokyo, 36, 99134, 1958. Acknowledgments: Raphaël Grandin, Benoît Taisne Download mogi.m at Matlab Central file exchange (Copyright 19972012 F. Beauducel / BSD License). 
The Sun [1969] model calculates analytical solution for surface deformation due to hydrostatic pressure inside a horizontal circular fracture (“pennyshaped”) in an elastic halfspace. The proposed Matlab script is a literal transcription of the Sun's paper equations. The equations are also vectorized for all input parameters. See help for syntax, and script comments for details. SUN69 Deformation from pennyshaped crack in elastic halfspace. [Ur,Uz] = SUN69(R,H,A,ΔV) or SUN69(R,H,A,ΔP,E,ν) computes radial and vertical displacements Ur and Uz on the free surface, due to a horizontal circular fracture formed in a semiinfinite elastic medium, with following variables: R: radial distance of observation, H: depth of the center of the source from the surface, A: radius of the source with the hydrostatic pressure, ΔV: volume of injected material, ΔP: change of the hydrostatic pressure in the crack. E: Young's modulus, ν: Poisson's ratio (default is 0.25 for isotropic medium). [Ur,Uz,B] = SUN69(...) returns also the maximum separation of fracture (vertical displacement) B. Equations from Sun [1969], with approximation H/A >> 1. If H/A > 2, error is about 2 to 3%; if H/A > 5, solution is almost perfect. Notes:  Equations are all vectorized, so variables R,H,A,ΔV or ΔP can be scalar or any of them vector or matrix, then outputs will have the same size.  Convention: Uz > 0 = UP, H is depth so in Z direction.  Units should be constistent, e.g.: R, H, A, Ur and Uz in m imply V in m³; optional E and ΔP in Pa, ν dimensionless. Example for a 3D plot of exagerated deformed surface: [x,y] = meshgrid(3:.1:3); [th,rho] = cart2pol(x,y); [ur,uz] = sun69(rho,1,.5,1e6,10e9,0.25); [ux,uy] = pol2cart(th,ur); ps = 5e4; surf(x+ux*ps,y+uy*ps,uz*ps), axis equal, light Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris, 2009. References: Sun, R. J. (1969). Theoretical size of hydraulically induced horizontal fractures and corresponding surface uplift in an idealized medium, J. Geophys. Res., 74, 59956011. Download sun69.m at Matlab Central file exchange (Copyright 2010 F. Beauducel / BSD License). 
The Okada [1985] model calculates analytical solution for surface deformation due to shear and tensile faults in an elastic halfspace. This model is widely used to simulate ground deformation produced by local perturbation like tectonic faults (earthquakes) or volcanic dykes (magmatic intrusion). Given rectangular fault geometry (length, width, depth, strike, dip) and 3component dislocation amplitude (rake, slip and open), it computes the displacements, tilts and strains at the freesurface. The proposed Matlab script is a literal transcription of the Okada's equations, except that it is transposed in a geographical referential (East, North, Up), where the fault is defined by a strike angle relative to the North, and dislocation parameters are given by: rake, slip and opening (instead of U1, U2, U3), following Aki & Richards [1980] definition. All coordinates and depth are relative to fault centroid. Lamé's constants λ and μ are replaced by Poisson's ratio ν (with a default value of 0.25 for isotropic medium), since the equations are independent of other elastic parameters. The equations are also vectorized for (x,y) coordinates and all input parameters except dip angle. To check the consistency of numerical calculations, run the script okada85_checklist.m, a transcription of table 2 cases 2, 3, and 4 checklist from [Okada, 1985] paper. See help for further details, syntax, example, and script comments for technical details. OKADA85 Surface deformation due to a finite rectangular source. [uE,uN,uZ,uZE,uZN,uNN,uNE,uEN,uEE] = OKADA85(... E,N,DEPTH,STRIKE,DIP,LENGTH,WIDTH,RAKE,SLIP,OPEN) computes displacements, tilts and strains at the surface of an elastic halfspace, due to a dislocation defined by RAKE, SLIP, and OPEN on a rectangular fault defined by orientation STRIKE and DIP, and size LENGTH and WIDTH. The fault centroid is located (0,0,DEPTH). E,N : coordinates of observation points in a geographic referential (East,North,Up) relative to fault centroid (units are described below) DEPTH : depth of the fault centroid (DEPTH > 0) STRIKE : fault trace direction (0 to 360° relative to North), defined so that the fault dips to the right side of the trace DIP : angle between the fault and a horizontal plane (0 to 90°) LENGTH : fault length in the STRIKE direction (LENGTH > 0) WIDTH : fault width in the DIP direction (WIDTH > 0) RAKE : direction the hanging wall moves during rupture, measured relative to the fault STRIKE (180 to 180°) SLIP : dislocation in RAKE direction (length unit) OPEN : dislocation in tensile component (same unit as SLIP) returns the following variables (same matrix size as E and N): uN,uE,uZ : displacements (unit of SLIP and OPEN) uZE,uZN : tilts (in rad * FACTOR) uNN,uNE,uEN,uEE : horizontal strains POSITIVE = COMPRESSION (unit of FACTOR) Length unit consistency: E, N, DEPTH, LENGTH, and WIDTH must have the same unit (e.g. km) which can be different from that of SLIP and OPEN (e.g. m) but with a possible FACTOR on tilt and strain results (in this case, an amplification of km/m = 1000). To have FACTOR = 1 (tilt in radians and correct strain unit), use the same length unit for all aforesaid variables. [...] = OKADA85(...,NU) specifies Poisson's ratio NU (default is 0.25 for an isotropic medium). Formulas and notations from Okada [1985] solution excepted for strain convention (here positive strain means compression), and for the fault parameters after Aki & Richards [1980], e.g.: DIP=90, RAKE=0 : left lateral (senestral) strike slip DIP=90, RAKE=180 : right lateral (dextral) strike slip DIP=70, RAKE=90 : reverse fault DIP=70, RAKE=90 : normal fault It is also possible to produce partial outputs, with following syntax: [uE,uN,uZ] = OKADA85(...) for displacements only; [uE,uN,uZ,uZE,uZN] = OKADA85(...) for displacements and tilts; [uE,uN,uZ,uNN,uNE,uEN,uEE] = OKADA85(...) for displacements and strains; [uZE,uZN] = OKADA85(...) for tilts only; [uZE,uZN,uNN,uNE,uEN,uEE] = OKADA85(...) for tilts and strains; [uNN,uNE,uEN,uEE] = OKADA85(...) for strains only. Note that vertical strain components can be obtained with following equations: uNZ = uZN; uEZ = uZE; uZZ = (uEE + uNN)*NU/(1NU); [...] = OKADA85(...,'plot') or OKADA85(...) without output argument produces a 3D figure with fault geometry and dislocation at scale (if all of the fault parameters are scalar). Equations are all vectorized excepted for argument DIP which must be a scalar (beacause of a singularity in Okada's equations); all other arguments can be scalar or matrix of the same size. Example: [E,N] = meshgrid(linspace(10,10,50)); [uE,uN,uZ] = okada85(E,N,2,30,70,5,3,45,1,1,'plot'); figure, surf(E,N,uN) considers a 5x3 fault at depth 2, N30° strike, 70° dip, and unit dislocation in all directions (reverse, senestral and open). Displacements are computed on a regular grid from 10 to 10, and North displacements are plotted as a surface. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 1997 Updated: 20140524 References: Aki K., and P. G. Richards, Quantitative seismology, Freemann & Co, New York, 1980. Okada Y., Surface deformation due to shear and tensile faults in a halfspace, Bull. Seismol. Soc. Am., 75:4, 11351154, 1985. Acknowledgments: Dmitry Nicolsky, Qian Yao, Halldor Geirsson. : new version on May 24, 2014 with bug corrections. Update your file! Download okada85.m at Matlab Central file exchange (Copyright 19972014 F. Beauducel / BSD License).

The Okubo [1992] model calculates analytical solution for gravity changes due to shear and tensile faults in an elastic halfspace (as Okada [1985] does for deformations). This model can be used to simulate gravity anomalies produced by local perturbation like tectonic faults (earthquakes) or volcanic dykes (magmatic intrusion). Given rectangular fault geometry (length, width, depth, strike, dip) and 3component dislocation amplitude (rake, slip and open), it computes the total gravity change (potential, cavity filling and freeair) at the freesurface. The proposed Matlab script is a literal transcription of the Okubo's equations, except for the fault geometry that it is transposed in a geographical referential (East, North, Up), where the fault is defined by a strike angle relative to the North, and dislocation parameters are given by: rake, slip and opening (instead of U1, U2, U3), following Aki & Richards [1980] definition. All coordinates and depth are relative to fault centroid. The equations are also vectorized for (x,y) coordinates and all parameters except dip angle. Type "doc okubo92" for help, syntax and example, and see script comments for details. OKUBO92 Surface gravity change due to a finite rectangular source. [dG,dH] = OKUBO92(E,N,DEPTH,STRIKE,DIP,LENGTH,WIDTH,RAKE,SLIP,OPEN,RHO,RHOP) computes total gravity change at the free surface of an elastic halfspace, due to a dislocation defined by RAKE, SLIP, and OPEN on a rectangular fault defined by orientation STRIKE and DIP, and size LENGTH and WIDTH. The fault centroid is located (0,0,DEPTH). E,N : coordinates of observation points in a geographic referential (East,North,Up) relative to fault centroid (units are described below) DEPTH : depth of the fault centroid (DEPTH > 0) STRIKE : fault trace direction (0 to 360° relative to North), defined so that the fault dips to the right side of the trace DIP : angle between the fault and a horizontal plane (0 to 90°) LENGTH : fault length in the STRIKE direction (LENGTH > 0) WIDTH : fault width in the DIP direction (WIDTH > 0) RAKE : direction the hanging wall moves during rupture, measured relative to the fault STRIKE (180 to 180°). SLIP : dislocation in RAKE direction (in m) OPEN : dislocation in tensile component (in m) RHO : density of the medium (in kg/m^3) RHOP : optional density of the cavityfilling matter (in kg/m^3); takes RHO value if not specified. returns the following variables (same matrix size as E and N): dG : total gravity change (in m/s^2) dH : elevation changes as given by Okada [1985] model (in m) Length unit consistency: E, N, DEPTH, LENGTH, and WIDTH must have the same unit (e.g. km) which can be different from that of SLIP and OPEN (e.g. m to to produce dG in m/s^2 and dH in m). [...] = OKUBO92(...,RHOP,BETA,NU) specifies gravity vertical gradient BETA (default is 0.309e5 s^2) and Poisson's ratio NU (default is 0.25 for an isotropic medium). Formulas and notations from Okubo [1992] solution excepted for the fault geometry parameters after Aki & Richards [1980], e.g.: DIP=90, RAKE=0 : left lateral (senestral) strike slip DIP=90, RAKE=180 : right lateral (dextral) strike slip DIP=70, RAKE=90 : reverse fault DIP=70, RAKE=90 : normal fault Equations are all vectorized excepted for argument DIP which must be a scalar; all other arguments can be scalar or matrix of the same size. Example: [E,N] = meshgrid(linspace(10,10,50)); dG = okubo92(E,N,6,90,90,10,10,0,5,0,2670); figure, pcolor(E,N,1e8*dG), shading interp hold on, [c,h]=contour(E,N,1e8*dG,50:10:50,'k'); clabel(c,h), hold off colorbar, polarmap reproduces the figure 4a in Okubo [1992] paper, considering a 10x10km fault at 6km depth, oriented N90°strike, 90°dip, and 5m senestral dislocation. Gravity changes are computed on a regular grid from 10 to 10km, and are plotted as a surface in 1e8 m/s^2 unit. POLARMAP is an author's colormap. See also OKADA85 author's function to compute complete displacements, tilt and strain at free surface. Author: Francois Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20120611 Updated: 20120612 References: Aki K., and P. G. Richards, Quantitative seismology, Freemann & Co, New York, 1980. Okubo S., Gravity and Potential Changes due to Shear and Tensile Faults in a HalfSpace, J. Geophys. Res., 97:B5, 71377144, 1992. Download okubo92.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 
The Standard for the Exchange of Earthquake Data (SEED) is an international standard format for the exchange of digital seismological data. SEED was designed for use by the earthquake research community, primarily for the exchange between institutions of unprocessed earth motion data. It is a format for digital data measured at one point in space and at equal intervals of time. The SEED format consists of Volume Control Headers, Abbreviation Control Headers, Station Control Headers, Time Span Control Headers and finally Data Records. In complement to “Dataless” SEED volumes, exists the “Dataonly” volume called MiniSEED (see http://www.iris.edu for further information). The purpose of these functions is to read and write miniSEED data files directly from Matlab, avoiding intermediate format conversion (like SAC or other formats for which many functions exist), having a full control on headers and formats. Note that IRIS proposes a complete Clibrary for miniSEED (libmseed) and there is a free Matlab toolbox using this library to read and write miniSEED files (see Stefan Mertl site), but it needs some compilation to be used (unsignificant for geeks, quite easy on UNIX/Linux, but no doubt more problematic for Windows users...). The present functions rdmseed.m and mkmseed.m are 100% Matlab code standalone scripts, slower, but fully portable and platform independent. rdmseed: reading miniSEED file Each data record is imported into a structure array, allowing to adress data blocks and header fields individually (useful for multichannel files), just as concatenating all data with a simple cat(1,X.d) function. Time stamps are also converted into Matlab datenum format. The function reads miniSEED "dataonly" using the two mostly used compression formats Steim1 and Steim2. General FDSN formats have also been implemented (ASCII, 16/24/32bit integers, IEEE floats and doubles), and GEOSCOPE multiplexed old formats (24bit, 16/3 or 16/4bit gain ranged). All these formats should work but some of them have not been tested using real data. I also partly coded Steim3 format but without a clear description and any file example... Since I never met any data file using this format, I don't know if it's really useful. The function detects also automatically big/littleendian coded files. Known Blockettes are 1000, 1001, 100, 500 and 2000. If there is no Blockette 1000 (which is mandatory in SEED format...), default 4096byte record length, bigendian and Steim1 compression are used. These values can be set using additional arguments. Using extra output argument, some analysis can be done on the data stream (detection of gaps and overlaps), and channel components are detected. Without any output arguments, or with an additionnal 'plot' input argument, the function plots the imported data in a new figure (works also in case of multichannel file). Steim1/2 compression decoding strategy has been deeply optimized for Matlab. The proposed method, as vectorized as possible, is about 30 times faster than a 'Clike' loops coding... which is still 10 times slower than the same Ccompiled program, but, well, this is the Matlab's other side of the coin! mkmseed: writing miniSEED file The function allows to export a data vector D to miniSEED file, giving origin date and time T0 and sampling rate FS (for monotonic data) or a time vector T. Header information is specified using the filename string with conventional naming "Network.Station.Location.Channel". Output file names will have appended ".Year.Day" and multiple file may be produced if data exceed a day. Data encoding format can be specified (16/32bit integers, IEEE float/double, Steim1/2; Geoscope 16/34). If not, it will depend on the class of variable D. Note that D will be converted to relevant class depending on encoding format. Be aware that Steim2 compression uses 30bit signed integer to code difference between 2 consecutive samples; equivalent to a 29bit coding. Steim1 uses 32bit signed integers to code differences, so a 31bit equivalent. This means that full range of signed 32bit integer data cannot be correctly encoded with these formats in some exceptional cases (may lead to data loss). Finally, Steim1/2 encoding needs time processing to achieve a maximum compression ratio of about 6 in file size, while other formats are very fast but result in large files. Binary file will be bigendian coded, with blockettes 1000, and default record length is 4096 bytes (this may be changed using input argument). Type "doc rdmseed" or "doc mkmseed" for detailed usage. RDMSEED Read miniSEED format file. X = RDMSEED(F) reads file F and returns a Mby1 structure X containing M blocks ("data records") of a miniSEED file with headers, blockettes, and data in dedicated fields, in particular, for each data block X(i): t: time vector (DATENUM format) d: data vector (double) BLOCKETTES: existing blockettes (substructures) Known blockettes are 100, 500, 1000, 1001 and 2000. Others will be ignored with a warning message. X = RDMSEED(F,ENCODINGFORMAT,WORDORDER,RECORDLENGTH), when file F does not include the Blockette 1000 (like Seismic Handler outputs), specifies:  ENCODINGFORMAT: FDSN code (see below); default is 10 = Steim1;  WORDORDER: 1 = bigendian (default), 0 = littleendian;  RECORDLENGTH: must be a power of 2, at least 256 (default is 4096). If the file contains Blockette 1000 (which is mandatory in the SEED convention...), these 3 arguments are ignored. X = RDMSEED without input argument opens user interface to select the file from disk. [X,I] = RDMSEED(...) returns a Nby1 structure I with N the detected number of different channels, and the following fields: ChannelFullName: channel name, XBlockIndex: channel's vector index into X, ClockDrift: vector of time interval errors, in seconds, between each data block (relative to sampling period). This can be compared to "Max Clock Drift" value of a Blockette 52. = 0 in perfect case < 0 tends to overlapping > 0 tends to gapping OverlapBlockIndex: index of blocks (into X) having a significant overlap with previous block (less than 0.5 sampling period). OverlapTime: time vector of overlapped blocks (DATENUM format). GapBlockIndex: index of blocks (into X) having a significant gap with next block (more than 0.5 sampling period). GapTime: time vector of gapped blocks (DATENUM format). RDMSEED(...) without output arguments plots the imported signal by concatenating all the data records, in one single plot if single channel is detected, or subplots for multichannels file. Gaps are shown with red stars, overlaps with green circles. [...] = RDMSEED(F,...,'be') forces bigendian reading (overwrites the automatic detection of endianness coding, which fails in some cases). [...] = RDMSEED(F,...,'plot') forces the plot with output arguments. [...] = RDMSEED(F,...,'v') uses verbose mode (displays additional information and warnings when necessary). Use 'vv' for extras, 'vvv' for debuging. Some instructions for usage of the returned structure:  to get concatenated time and data vectors from a singlechannel file: X = rdmseed(f,'plot'); t = cat(1,X.t); d = cat(1,X.d);  to get the list of channels in a multichannel file: [X,I] = rdmseed(f); cat(1,I.ChannelFullName)  to extract the station component i from a multichannel file: [X,I] = rdmseed(f); k = I(i).XBlockIndex; plot(cat(1,X(k).t),cat(1,X(k).d)) datetick('x') title(I(i).ChannelFullName) Known encoding formats are the following FDSN codes: 0: ASCII 1: 16bit integer 2: 24bit integer (untested) 3: 32bit integer 4: IEEE float32 5: IEEE float64 10: Steim1 11: Steim2 12: GEOSCOPE 24bit (untested) 13: GEOSCOPE 16/3bit gain ranged 14: GEOSCOPE 16/4bit gain ranged 19: Steim3 (alpha and untested) See also MKMSEED to export data in miniSEED format. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20100917 Updated: 20140507 Acknowledgments: Ljupco Jordanovski, JeanMarie Saurel, Mohamed Boubacar, Jonathan Berger, Shahid Ullah, Wayne Crawford, Constanza Pardo. References: IRIS (2010), SEED Reference Manual: SEED Format Version 2.4, May 2010, IFDSN/IRIS/USGS, http://www.iris.edu Trabant C. (2010), libmseed: the MiniSEED library, IRIS DMC. Steim J.M. (1994), 'Steim' Compression, Quanterra Inc. MKMSEED Write data in miniSEED file format. MKMSEED(FILENAME,D,T0,FS) writes miniSEED file FILENAME from stritl monotonic data vector D, time origin T0 (a scalar in Matlab datenum compatible format) and sampling rate FS (in Hz). Encoding format will depend on D variable class (see below). MKMSEED(FILENAME,D,T,FS) where T is a time vector of the same length as data vector D, will create data records of monotonic blocks of samples, splitting each time the sampling frequency FS is not equal to time difference between two successive samples (with a 50% tolerency). Network, Station, Channel and Location codes will be extracted from FILENAME which must respect the format "NN.SSSSS.LC.CCC" where: NN = Network Code (2 characters max, see FDSN list) SSSSS = Station identifier (5 char max) LC = Location Code (2 char max) CCC = Channel identifier (3 char max) Final filename will have appended string ".YYYY.DDD" corresponding to year and ordinal day of origin time (from T0 value). Multiple files may be created if time span of data exceeds day limit. MKMSEED(...,EF,RL) specifies encoding format EF and data record length RL (in bytes). RL must be a power of 2 greater or equal to 256. Data encoding format EF must be one of the following FDSN codes: 1 = 16bit integer (default for class 2bit, 8bit, int16) 3 = 32bit integer (default for class uint16, int32) 4 = IEEE float32 (default for class single) 5 = IEEE float64 (default for all other class) 10 = Steim1 compression (D will be converted to int32) 11 = Steim2 compression (D will be converted to int32) 13 = Geoscope 163 gain ranged (D will be converted to double) 14 = Geoscope 164 gain ranged (D will be converted to double) MKMSEED(FILENAME,D) uses default value for T0 = now (present date and time), and FS = 100 Hz. File(s) will be coded bigendian, flags set to zero, blockette 1000, default record length is 4096 bytes. Outputs have been tested with PQL II software from IRIS PASSCAL (v2010.268). See also RDMSEED function for miniSEED file reading. Acknowledgments: Florent Brenguier, Julien Vergoz, Constanza Pardo. References: IRIS (2010), SEED Reference Manual: SEED Format Version 2.4, May 2010, IFDSN/IRIS/USGS, http://www.iris.edu IRIS (2010), PQL II Quick Look trace viewing application, PASSCAL Instrument Center, http://www.passcal.nmt.edu/ Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20111019 Updated: 20140507 : new version of rdmseed on April 21, 2012 with bug corrections (littleendian coding). Update your file! Download rdmseed.m and mkmseed.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). 
RDSAC reads a seismic data file encoded in the IRIS/SAC binary format, and returns a time vector, a data vector and all header variables in a structure (field names correspond to exact IRIS variable names). Type "doc rdsac" for detailed usage. RDSAC Read SAC data file. X=RDSAC(FILE) reads the Seismic Analysis Code (SAC) FILE and returns a structure X containing the following fields: t: time vector (DATENUM format) d: data vector (double) HEADER: header substructure (as defined in the IRIS/SAC format). RDSAC without input argument will open a file browser window. X=RDSAC(...,'plot') or RDSAC(...) without output argument will plot the data in a new figure. Reference: http://www.iris.edu/files/sacmanual/ Author: F. Beauducel <beauducel@ipgp.fr> Created: 20140401 Updated: 20140425 Download rdsac.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). 
This function plots regular grids of elevation X,Y,Z in a more efficient manner than SURFL Matlab's function, because it recomputes lighting and displays result as shaded color flat RGB image. It uses also medianstyle filter to deal with min/max values of elevations and gradient, and proposes two specific colormaps "landcolor" and "seacolor". Color mapping and lighting parameters can be changed from default values. In addition, several options are available: 'cartesian' to add decimal axis, 'latlon' to add geographical axis (GMTlike), 'legend' for an automatic scaling legend, 'lake' for automatic flat area colorfilling and 'interp' to fill the novalue gaps... This may be useful to produce highquality and moderatesize Postscript image adapted for publication. Figure examples:
See "doc dem" for syntax, examples and help. See also the READHGT function from same author. DEM Shaded relief image plot DEM(X,Y,Z) plots the Digital Elevation Model defined by X and Y coordinate vectors and elevation matrix Z, as a lighted image using specific "landcolor" and "seacolor" colormaps. DEM uses IMAGESC function which is much faster than SURFL when dealing with large highresolution DEM. It produces also highquality and moderatesize Postscript image adapted for publication. [H,I] = DEM(...); returns graphic handle H and illuminated image as I, an MxNx3 matrix (if Z is MxN and DECIM is 1). DEM(X,Y,Z,'Param1',Value1,'Param2',Value2,...) specifies options or parameter/value couple (case insensitive):  Lighting options  'Azimuth',A Light azimuth in degrees clockwise relative to North. Default is A = 45 for a natural northwestern illumination. 'Contrast',C Light contrast, as the exponent of the gradient value: C = 1 for linear contrast (default), C = 0 to remove lighting, C = 0.5 for moderate lighting, C = 2 or more for strong contrast. 'LCut',LC Lighting scale saturation cut with a medianstyle filter in % of elements, such as LC% of maximum gradient values are ignored: LC = 0.2 is default, LC = 0 for full scale gradient. 'km' Stands that X and Y coordinates are in km instead of m (default). This allows correct lighting. Ignored if LATLON option is used.  Elevation colorscale options  'ZLim',[ZMIN,ZMAX] Fixes min and max elevation values for colormap. Use NaN to keep real min and/or max data values. 'ZCut',ZC Medianstyle filter to cut extremes values of Z (in % of elements), such that ZC% of most min/max elevation values are ignored in the colormap application: ZC = 0.5 is default, ZC = 0 for full scale.  "No Value" elevation options  'NoValue',NOVALUE Defines the values that will be replaced by NaN. Note that values equal to minimum of Z class are automatically detected as NaN (e.g., 32768 for int16 class). 'NaNColor',[R,G,B] Sets the RGB color for NaN/NoValue pixels (default is a dark gray). Note that your must specify a valid 3scalar vector (between 0 and 1); color characters like 'w' or 'k' are not allowed, use [1,1,1] or [0,0,0] instead. 'Interp' Interpolates linearly all NaN values (fills the gaps using linear triangulation), using an optimized algorithm.  Colormap options  'LandColor',LMAP Uses LMAP colormap instead of default (landcolor, if exists or jet) for Z > 0 elevations. 'SeaColor',SMAP Sets the colormap used for Z ≤ 0 elevations. Default is seacolor (if exists) or single color [0.7,0.9,1] (a light cyan) to simulate sea color. 'ColorMap',CMAP Uses CMAP colormap for full range of elevations, instead of default land/sea. This option overwrites LANDCOLOR/SEACOLOR options. 'Lake' Detects automaticaly flat areas different from sea level (nonzero elevations) and colors them as lake surfaces. 'Watermark',N Makes the whole image lighter by a factor of N.  Basemap and scale options  'Cartesian' Plots classic basemapstyle axis, considering coordinates X and Y as cartesian in meters. Use parameter "km' for X/Y in km. 'LatLon' Plots geographic basemapstyle axis in deg/min/sec, considering coordinates X as longitude and Y as latitude. Axis aspect ratio will be adjusted to approximatively preserve distances (this is not a real projection!). This overwrites ZRatio option. 'Legend' Adds legends to the right of graph: elevation scale (colorbar) and a distance scale (in km). 'FontSize',FS Font size for X and Y tick labels. Default is FS = 10. 'BorderWidth',BW Border width of the basemap axis, in % of axis height. Must be used together with CARTESIAN or LATLON options. Default is BW = 1%.  Decimation options  For optimization purpose, DEM will automatically decimate data to limit to a total of 1500x1500 pixels images. To avoid this, use following options, but be aware that large grids may require huge computer ressources or induce disk swap or memory errors. 'Decim',N Decimates matrix Z at 1/N times of the original sampling. 'NoDecim' Forces full resolution of Z, no decimation.  Informations  Colormaps are Mx3 RGB matrix so it is easy to modify saturation (CMAP.^N), set darker (CMAP/N), lighter (1  1/N + CMAP/N), inverse it (flipud(CMAP)), etc... To get free worldwide topographic data (SRTM), see READHGT function. For backward compatibility, the former syntax is still accepted: DEM(X,Y,Z,OPT,CMAP,NOVALUE,SEACOLOR) where OPT = [A,C,LC,ZMIN,ZMAX,ZC], also option aliases DEC, DMS and SCALE, but there is no argument checking. Please prefer the param/value syntax. Author: François Beauducel <beauducel@ipgp.fr> Created: 20070517 Updated: 20130310 Download dem.m, landcolor.m and seacolor (needed) at Matlab Central file exchange (Copyright 2013 F. Beauducel / BSD License). 
This function imports .HGT "height" binary data files from NASA SRTM global digital elevation model of Earth land, corresponding to 1x1 degree tiles of 3arc seconds resolution (SRTM3, around 90 m) and 1arc second (SRTM1, around 30 m) for USA territory, and returns coordinates vectors latitude and longitude, and a matrix of elevation values. The function includes an automatic download of data from the USGS SRTM webserver, so indicating latitude and longitude is sufficient to get the data and instant map plot anywhere in the World. Also some basic options as 'merge' to concatenate tiles, 'interp' for gaps (novalue) linear interpolation, 'crop' for rectangle selection inside tile(s). READHGT Import/download NASA SRTM data files (.HGT). X=READHGT(FILENAME) reads HGT "height" SRTM data file FILENAME and returns a structure X containing following fields lat: coordinate vector of latitudes (decimal degree) lon: coordinate vector of longitudes (decimal degree) z: matrix of elevations (meters, INT16 class) FILENAME must be in the form "[NS]yy[EW]xxx.hgt[.zip]" as downloaded from SRTM data servers. X=READHGT(LAT,LON) attemps to download the file *.hgt.zip corresponding to LAT and LON (in decimal degrees) coordinates (lowerleft corner) from the USGS data server (needs an Internet connection and a companion file "readhgt_srtm_index.txt"). Downloaded filename(s) will be given in an additional output structure field X.hgt. LAT and/or LON can be vectors: in that case, tiles corresponding to all possible combinations of LAT and LON values will be downloaded, and optional output structure X will have as much elements as tiles. READGHT(...) without output argument or X=READHGT(...,'plot') plots the tile(s). For better plot results, it is recommended to install DEM personal function available at author's Matlab page. READHGT(LAT,LON,...,'merge'), in case of adjoining values of LAT and LON, will concatenate tiles to produce a single one. ATTENTION: tiles assembling may require huge computer ressources and cause disk swaping or memory error. READHGT(LAT,LON,...,'interp') linearly interpolates missing data. READHGT(...,'crop',[LAT1,lAT2,LON1,LON2]) crops the map using latitude/longitude limits. READHGT(LAT,LON,...,'crop'), without limits argument vector, crops the resulting map around existing land (reduces any sea or novalue areas at the borders). READHGT(LAT,LON,...,'srtm3') forces SRTM3 download (by default, SRTM1 tile is downloaded if exists). Usefull for USA neighborhood. READHGT(LAT,LON,OUTDIR) specifies output directory OUTDIR to write downloaded files. READHGT(LAT,LON,OUTDIR,URL) specifies the URL address to find HGT files (default is USGS). Examples:  to plot a map of the Paris region, France (single tile): readhgt(48,2)  to plot a map of Flores volcanic island, Indonesia (5 tiles): readhgt(9,119:123,'merge')  to download SRTM1 data of Cascade Range (27 individual tiles): X=readhgt(40:48,123:121); Informations:  each file corresponds to a tile of 1x1 degree of a square grid 1201x1201 of elevation values (SRTM3 = 3 arcseconds), and for USA territory at higher resolution 3601x3601 grid (SRTM1 = 1 arcsecond)  elevations are of class INT16: sea level values are 0, unknown values equal 32768 (there is no NaN for INT class), use 'interp' option to fill the gaps  note that borders are included in each tile, so to concatenate tiles you must remove one row/column in the corresponding direction (this is made automatically with the 'merge' option)  downloaded file is written in the current directory or optional OUTDIR directory, and it remains there  NASA Shuttle Radar Topography Mission [February 11 to 22, 2000] produced a nearglobal covering on Earth land, but still limited to latitudes from 60S to 60N. Offshore tiles will be output as flat 0 value grid  if you look for other global topographic data, take a look to ASTER GDEM, worldwide 1 arcsecond resolution (from 83S to 83N): http://gdem.ersdac.jspacesystems.or.jp (free registration required) Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris References: http://dds.cr.usgs.gov/srtm/version2_1 Acknowledgments: Yves Gaudemer Created: 20120422 Updated: 20130117 Download readhgt.m and readhgt_srtm_index.txt (needed) at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 
UTM2LL converts Universal Transverse Mercator (UTM) East/North coordinates to latitude/longitude. LL2UTM converts latitude/longitude coordinates to UTM. Both functions are using precise formula (millimeter precision), possible userdefined datum (WGS84 is the default), and are all vectorized (no loop in the code). It means that huge matrix of points, like an entire DEM grid, can be converted very fast. Example (needs readhgt.m author's function): X = readhgt(36:38,12:15,'merge','crop',[36.5,38.5,12.2,16],'plot'); [lon,lat] = meshgrid(X.lon,X.lat); [x,y,zone] = ll2utm(lat,lon); z = double(X.z); z(z==32768  z<0) = NaN; figure pcolor(x,y,z); shading flat; hold on contour(x,y,z,[0,0],'w') hold off; axis equal; axis tight xlabel('East (m)'); ylabel('North (m)') title(sprintf('Sicily  UTM zone %d WGS84',zone)) loads SRTM full resolution DEM of Sicily in lat/lon (a 2400x4500 grid), converts it to UTM and plots the result. To make a regular UTM grid, you may interpolate x and y with griddata function. See "doc ll2utm" and "doc utm2ll" for syntax and help. LL2UTM Lat/Lon to UTM coordinates precise conversion. [X,Y]=LL2UTM2(LAT,LON) converts coordinates LAT,LON (in degrees) to UTM X and Y (in meters). Default datum is WGS84. LAT and LON can be scalars, vectors or matrix. Outputs X and Y will have the same size as inputs. LL2UTM(LAT,LON,DATUM) uses specific DATUM for conversion. DATUM can be a string in the following list: 'wgs84': World Geodetic System 1984 (default) 'nad27': North American Datum 1927 'clk66': Clarke 1866 'nad83': North American Datum 1983 'grs80': Geodetic Reference System 1980 'int24': International 1924 / Hayford 1909 or DATUM can be a 2element vector [A,F] where A is semimajor axis (in meters) and F is flattening of the userdefined ellipsoid. [X,Y,ZONE]=LL2UTM(...) returns also computed UTM ZONE (negative value stands for southern hemisphere points). Notice:  LL2UTM does not perform crossdatum conversion.  precision is near a millimeter. Reference: I.G.N., Projection cartographique Mercator Transverse: Algorithmes, Notes Techniques NT/G 76, janvier 1995. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20031202 Updated: 20140217 UTM2LL UTM to Lat/Lon coordinates precise conversion. [LAT,LON]=UTM2LL(X,Y,ZONE) converts UTM coordinates X,Y (in meters) defined in the UTM ZONE (integer) to latitude LAT and longitude LON (in degrees). Default datum is WGS84. X and Y can be scalars, vectors or matrix. Outputs LAT and LON will have the same size as inputs. For southern hemisphere points, use negative zone ZONE. UTM2LL(X,Y,ZONE,DATUM) uses specific DATUM for conversion. DATUM can be a string in the following list: 'wgs84': World Geodetic System 1984 (default) 'nad27': North American Datum 1927 'clk66': Clarke 1866 'nad83': North American Datum 1983 'grs80': Geodetic Reference System 1980 'int24': International 1924 / Hayford 1909 or DATUM can be a 2element vector [A,F] where A is semimajor axis (in meters) and F is flattening of the userdefined ellipsoid. Notice:  UTM2LL does not perform crossdatum conversion.  precision is near a millimeter. Reference: I.G.N., Projection cartographique Mercator Transverse: Algorithmes, Notes Techniques NT/G 76, janvier 1995. Author: Francois Beauducel, <beauducel@ipgp.fr> Created: 20010823 Updated: 20140217 Download ll2utm.m and utm2ll.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). 

This program computes a map of visibility from a selected point on a topography. It has been written to help the search for radio relay best location. Because it considers only direct line of sight, it gives a good estimation for possible radio link for short distances only (less than 10 km), neglecting curvature of the Earth, Fresnel zone and atmospheric refraction on radio waves propagation. The program computes the relative elevation angle of the mask for each point (the angle is null or negative if the point is visible). The function needs a digital elevation model Z and associated (X,Y) vectors or matrices of coordinates (same unit as Z), position of the point (X0,Y0), the antenna height H0 (for instance 4 m), and the hypothetic antenna height Ha on each topography points (for instance 3 m). When no output argument is given, the function plots a map of the results (color map of mask angles, and blank for visible points, see example screenshot). The script is not fully optimized because it makes a global loop on the matrix elements to compute each profile of topography... (I didn't find (yet) a way to fully vectorize the problem), so it takes some time to compute, depending on the number of element of Z. However, I found a faster algorithm (about 2 times faster), giving approximate result but usefull to process a first map. See help for syntax, and script comments for details. RADIOCOVER Radio link coverage map. RADIOCOVER(X,Y,Z,X0,Y0,H0,Ha) computes the coverage map of possible direct linear radio link from the point (X0,Y0) with antenna height H0, using digital terrain model defined by coordinate vectors X and Y, and elevation matrix Z, and hypothetic antenna height Ha, then plots a color map of the relative elevation mask angle (in degrees) with blank areas where there is no mask (visible), together with a contour map of the topography. X and Y can be vectors with length(X) = n and length(Y) = m where [m,n] = size(Z), or matrices of the same size as Z (as from MESHGRID). RADIOCOVER(...,METHOD) specifies alternate methods. The default is nearest neighbor interpolation. Available methods are: 'nearest'  nearest neighbor interpolation (default) 'linear'  linear interpolation (smoother result) 'fast'  approximate algorithm (about 2 times faster) M = RADIOCOVER(...); returns a matrix of relative elevation mask angle (in degrees, same size as Z), without producing graphic. Visible points have null or negative values. The model assumes linear propagation of radio waves (direct line of sight between the two antennas), and neglects curvature of the Earth, Fresnel zone, and atmospheric refraction. Example: [x,y,z]=peaks(100); [fx,fy]=gradient(z); z=sqrt(fx.^2+fy.^2); surf(x,y,z), shading flat, light, view(24,74) radiocover(x,y,z,0.84,0.27,.05,.05) Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20030110 Modified: 20100108 Download radiocover.m at Matlab Central file exchange (Copyright 20032009 F. Beauducel / BSD License). 

The function GREATCIRCLE computes the shortest path along the great circle ("as the crow flies") between two points defined by their geographic coordinates (latitude and longitude). With one output argument it returns distance or vector of distances, with two or more output arguments it returns path coordinates and optional vector of distances and bearing angles along the path. The function LOXODROME computes the path with a constant bearing, crossing all meridians of longitude at the same angle. It returns also a vector of distances and the bearing angle. Loxodrome path (also known as "rhumb line") is longer than great circle one, but still used in navigation as it is easier to follow with a compass. Type 'doc greatcircle' or 'doc loxodrome' for syntax, help and examples. GREATCIRCLE "As the crow flies" path, distance and bearing. GREATCIRCLE(LAT1,LON1,LAT2,LON2) returns the shortest distance (in km) along the great circle between two points defined by spherical coordinates latitude and longitude (in decimal degrees). If input arguments are vectors or matrix, it returns a vector/matrix of distances of the same size. [LAT,LON,DISTANCE,BEARING]=GREATCIRCLE(LAT1,LON1,LAT2,LON2) where all input arguments are scalars, computes the way points of coordinates LAT and LON, the distances (in km), and the bearing angles (in degrees from North) along the path, with default 100 intermediate points. Note that last bearing angle is NaN. [...]=GREATCIRCLE(...,N) uses N intermediate points. Example: load topo contour(0:359,89:90,topo,[0,0],'k') [lat,lon,dis] = greatcircle(48.8,2.3,35.7,139.7); hold on, plot(lon,lat,'r','linewidth',2), hold off title(sprintf('Paris to Tokyo = %g km',dis(end))) See also LOXODROME. Author: Francois Beauducel <beauducel@ipgp.fr> Created: 20120726 Updated: 20121114 References: http://www.movabletype.co.uk/scripts/latlong.html http://en.wikipedia.org/wiki/Haversine_formula LOXODROME Rhumb line path and distance. [LAT,LON]=LOXODROME(LAT1,LON1,LAT2,LON2) computes the line between two points defined by their spherical coordinates latitude and longitude (in decimal degrees), crossing all meridians of longitude at the same angle, i.e. a path derived from an initial and constant bearing. LAT and LON contain vectors of coordinates of the way points. [LAT,LON,DISTANCE,BEARING]=LOXODROME(...) returns also vector of distances (in km) along the path way, and constant bearing angle (in degrees from North, a scalar). [...]=LOXODROME(...,N) uses N intermediate points (default is 100). Example: load topo contour(0:359,89:90,topo,[0,0],'k') [lat,lon,dis,bear] = loxodrome(48.8,2.3,35.7,139.7); hold on, plot(lon,lat,'r','linewidth',2), hold off title(sprintf('Paris to Tokyo = %g km  bear = %g N',dis(end),bear)) See also GREATCIRCLE. Author: Francois Beauducel <beauducel@ipgp.fr> Created: 20121030 Updated: 20121108 References: http://en.wikipedia.org/wiki/Rhumb_line Acknowledgments: Jorge David Taramona Perea, Jacques Vernin. Download greatcircle.m and loxodrome.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 

This little function draws a compass rose with optional Cardinal directions with basic properties:
COMPROSE Compass rose plot COMPROSE(X,Y,N,W,AZ) adds a compass rose on current axis located at position X,Y with N points (N is 1, 4, 8 or 16), width W (radius) and North pointing to azimuth AZ (in degree, AZ = 0 means an arrow pointing to positive Yaxis direction, rotating clockwise). COMPROSE(...,FS) adds Cardinal directions with font size FS. COMPROSE(...,'param1',value1,'param2',value2,...) specifies any additionnal properties of the Patch using standard parameter/value pairs. Note that 'FaceColor' concerns only the default blackfilled parts of the drawing. H=COMPROSE(...) returns a Nx3 matrix of object handles: first column addresses solid filled patches, second column for white patches, third column for Cardinal direction text. Examples: comprose(0,0,8,.5,10) comprose(2,1,1,2,0,20,'LineWidth',2,'FaceColor',.5*[1,1,1]) h = comprose(1,2.5,16,1,10); set(h(:,1),'FaceColor',[0,.5,.5],'EdgeColor',.5*[1,1,1]) set(h(:,2),'FaceColor',.9*[1,1,1]) See also PATCH. Author: Francois Beauducel <beauducel@ipgp.fr> Created: 20120624 Download comprose.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 

This little function could be an alternative to other existing arrow plot functions, since it has very simple, vectorized and effective coding. In fact, arrows are all plotted in a single patch command! The function also accepts any standard patch parameter/value pair like 'FaceColor', 'EdgeColor', 'LineWidth', etc. It can be used to plot a personalized arrow at coordinates X,Y with length L and azimuth AZ, or any of these arguments can be vector or matrix, like QUIVER function, to plot multiple arrows. Arrow's aspect ratio, head and segment line shapes are configurable with 4 parameters: head width, head length, head inside length and segment line width, all normalized to arrow's length. Type 'doc arrows' for syntax, help and examples. ARROWS Generalised 2D arrows plot ARROWS(X,Y,L,AZ) draws an arrow on current axis from position X,Y with length L, oriented with azimuth AZ (in degree, AZ = 0 means an arrow pointing to positive Yaxis direction, rotating clockwise). X and Y can be scalars or matrix. In the last case, any or both L and AZ can be scalars or matrix of the same size as X and Y. ARROWS(...,SHAPE) uses relative ratios SHAPE = [HEADW,HEADL,HEADI,LINEW] to adjust head width HEADW, head length HEADL, head inside length HEADI, and segment line width LINEW for an arrow length of 1 (default is SHAPE = [1,0.25,0.25,0.5]). ARROWS(...,'param1',value1,'param2',value2,...) specifies any additionnal properties of the Patch using standard parameter/value pairs, like 'FaceColor','EdgeColor','LineWidth', ... H=ARROWS(...) returns graphic's handle of patches. Examples: arrows(0,0,1,45,'FaceColor','none','LineWidth',3) arrows(1,0,1,0,[.2,.4,.2,.02]) [xx,yy] = meshgrid(1:10); arrows(xx,yy,rand(size(xx)),360*rand(size(xx))) Notes: Arrow shape supposes an equal aspect ratio (axis equal). Equivalent of quiver(x,y,u,v,0) command is obtained with: [th,l] = cart2pol(u,v); arrows(x,y,l,90  th*180/pi) See also PATCH, QUIVER. Author: Francois Beauducel <beauducel@ipgp.fr> Created: 19950203 Updated: 20120630 Download arrows.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 
This function gives the parameters for about 40 main wave components of Earth tides due to Sun and Moon (annual, semiannual, monthly, fortnightly, diurnal, semidiurnal, ...) using the harmonic development of Arthur Thomas Doodson [18901968]. For a given wave, it returns in a structure the 6 Doodson's arguments, wave Darwin's symbol, wave description and the wave period (in days). >> X=doodson('M2') X = symbol: {'M2'} name: {'Principal lunar semidiurnal'} doodson: [2 0 0 0 0 0] period: 0.5175Type doodson without argument to see the list of available waves. DOODSON Doodson tidal wave components X = DOODSON(WAVE) returns a structure containing, for the tidal wave symbol WAVE, the following fields: doodson: vector of the 6 Doodson arguments symbol: wave symbol string (Darwin's notation) name: wave long name period: wave period (in days) X = DOODSON returns all known waves parameters. Type DOODSON without any input/output argument displays a table of available waves. References: Agnew D.C. (2007), Earth Tides, in « Treatise on Geophysics: Geodesy », T. A. Herring Ed., Elsevier, New York. Author: Francois Beauducel <beauducel@ipgp.fr> Created: 20140522 Updated: 20140524 Download doodson.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). 
This function reads data files from Fluke TI32 infrared thermal camera (.IS2 proprietary format). It returns data images (raw and calibrated temperatures, and visible RGB image) and some of the header information in a single structure. Without output arguments, it plots blended infrared and visible images and a histogram of temperatures. The function has been tested for TI32 camera original files only. Presently it does not work with files edited with the Fluke software, neither for other cameras. Also the function needs a separated calibration file since I was not able to find the original tables in the binary file (but I'm quite sure they exist...). Any help to improve that is welcome. READIS2 Import data files from Fluke infrared thermal camera (IS2 format). X=READIS2(FILENAME) reads IS2 radiometric file FILENAME and returns a structure X containing fields with data and headers: IR: infrared sensor values (240x320 uint16 matrix) TB: brillancy temperatures, in °C (240x320 double matrix) T: body temperatures, in °C (240x320 double matrix) VI: visible image (480x640x3 uint8 RGB) VI2: cropped and interpolated visible image adjusted on IR image pixels (240x320x3 uint8 RGB) Emissivity: emissivity value (0 to 1.00) Transmission: transmission value (0 to 1.00) BackgroundTemp: background temperature (°C) Date: [YEAR MONTH DAY HOUR MINUTE SECOND] X=READIS2(FILENAME,MAP) plots temperature image using MAP color scale over visible image. MAP is a valid [n x 3] RGB colormap, for instance JET or GRAY. Use MAP=[] to specify other parameters without plot. X=READIS2(FILENAME,MAP,PARAM) applies given parameters to overwrite values from camera. PARAM = [DT,EM,TR,BT] (use NaN to keep original values): DT = time difference (in days) applied to image timestamp EM = emissivity TR = transmission BT = background temperature (°C) X=READIS2(FILENAME,MAP,PARAM,CAL) applies a polynomial function to calibrate brillancy temperatures, as follow: TB_cal = polyval(CAL,TB_raw). This will affect also body temperatures. Example: CAL = [1.08,2.72] is for my own Ti32 camera calibration result for 0100°C range laboratory experiment. Notes:  the function has been tested only with original IS2 files from a Fluke TI32 camera.  to produce temperatures from IR sensor values, the function needs a calibration file 'readis2_ti32.dat' located in current directory or Matlab path (see corresponding header info).  body temperatures Tr (output X.T) and brillancy temperatures Tb (output X.IR) are related with the formula: Tb = ε*τ*Tr + (2ετ)*Te where ε is emissivity, τ is transmission, and Te is background (reflected) temperature. So Tb and Tr are equal for ε = τ = 1. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20110724 Updated: 20120401 Acknowledgments: Magnus Andersen (read_is2.m), Pierre Agrinier, Valérie Clouard, Damien Gaudin, David Krammer. Download readis2.m and readis2_ti32.dat (needed) at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 
When representing "polarized" data for which zero and sign are meaningful (like for instance, a gravity anomalies map), it is useful to have a colormap where white stands for zero value, and greater magnitudes are plotted as colors with increasing saturation. This little function POLARMAP has two different usages:
POLARMAP Polarized color map POLARMAP applies a "polarized" redwhiteblue colormap to current figure, and adjusts the color axis limits to be centered to zero. POLARMAP(M) fixes the number of colors to M (default is 64). POLARMAP(MAP) applies linear shading to white to the center of colormap MAP which can be any of existing colormaps (an Mx3 matrix of RGB). POLARMAP(MAP,C) uses exponent C to modify the shading contrast. Default is C = 1 for linear shading. Use C = 2 to strengthen the shading, or C = 0.5 to attenuate it. C=POLARMAP(...) returns an Mby3 matrix containing the colormap, that can be used with COLORMAP function like other colormaps. Examples: pcolor(peaks), shading interp polarmap, colorbar then try the following polarmap(jet,0.5) Note the polar shading has no real interest with colormaps that include white color as one of the extremes (like GRAY, BONE, HOT, ...). See also JET, HSV, COPPER, SPRING, SUMMER, WINTER, COOL, COLORMAP, RGBPLOT. Author: François Beauducel <beauducel@ipgp.fr> Created: 20111026 Updated: 20120609 Download polarmap.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 
This function is a simple digital filter with one's coefficients, nearly equivalent to moving average but instead of mean values, it computes the sum. This is particularly usefull for rainfall plots for which it is, in my opinion, the only good way to show the data. The figure attempts to demonstrate it. Let us have some rainfall data at 10minute sampling interval, and simply want to plot the 1hour rainfall... If you use classical histogram representation, it leads to artifacts that depends on a priori time windows offset you chose (see figure):
MOVSUM Moving sum plot. Y = MOVSUM(X,N) returns the "moving sum" of X on N consecutive values. Y has the same size as X, and each element of Y is the sum of the N previous values of X. If X is a regularly spaced data vector, then Y corresponds to a digital filter with N bcoefficient equal to one: y(n) = x(n) + x(n1) + ... + x(1) This is nearly equivalent to a moving average, but instead of mean it is the sum. If X is a matrix, MOVSUM works down the columns. A typical use is for rainfall data vectors. For instance, if X contains rainfall values at 1minute interval, MOVSUM(X,60) will be the hourly rainfall, expressed as a continous vector of the same size as X. Compared to classical histogram representation, this method avoids some artifacts due to windowing using a priori time intervals. MOVSUM(...) without output arguments produces a plot of the results. Author: François Beauducel <beauducel@ipgp.fr> Created: 2005 Updated: 20101010 Download movsum.m (Copyright 2010 F. Beauducel / BSD License). 
This little function rounds a number (or the elements of a vector or matrix) towards the nearest number with N significant digits. This is a useful complement to Matlab's ROUND, ROUND10 and ROUNDN (Mapping toolbox), especially when dealing with data with a large variety of order of magnitudes. ROUNDSD Round with fixed significant digits ROUNDSD(X,N) rounds the elements of X towards the nearest number with N significant digits. ROUNDSD(X,N,METHOD) uses following methods for rounding: 'round'  nearest (default) 'floor'  towards minus infinity 'ceil'  towards infinity 'fix'  towards zero Examples: roundsd(0.012345,3) returns 0.0123 roundsd(12345,2) returns 12000 roundsd(12.345,4,'ceil') returns 12.35 See also Matlab's functions ROUND, ROUND10, FLOOR, CEIL, FIX, and ROUNDN (Mapping Toolbox). Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Acknowledgments: Edward Zechmann, Daniel Armyr, Yuri Kotliarov Created: 20090116 Updated: 20150403 Download roundsd.m at Matlab Central file exchange (Copyright 2015 F. Beauducel / BSD License). NEWS: This file was selected as MATLAB Central Pick of the Week, and integrated as an option of the core function ROUND in the new Matlab 2014b ! 
A little function that I imagined many years ago: a generalized "medianstyle" filter, an elegant way to generalize MIN, MAX, MEDIAN and extreme filtering functions, reduced to a single parameter. Let's sorting all values of a vector X in ascending order, thus consider the normalized rank position N of elements:
A typical example of use is minmax(X,[0.01 0.99]) which returns minimum and maximum values of X(:) but excluding the 1% amount of extreme values. This is particularily efficient for automatic scaling of noisy data (see the screenshot example), compared to the use of MEAN and STD functions which is biased by any highmagnitude values in X. Default behavior of MINMAX is to return a vector of minimum and maximum values. So minmax(X) is an abreviation of minmax(X,[0 1]), and equivalent of [min(X(:)) max(X(:))]. Compared to the builtin functions MIN and MAX, it does not return a vector from matrix but returns a scalar by considering all elements of matrix X. Also, it removes any NaN values (returns NaN only if all X values are NaN). Type 'doc minmax' for syntax, help and other examples. MINMAX Generalized median filter. Y=MINMAX(X) returns a 2 element vector Y = [min(X) max(X)]. Minimum and maximum values of X(:) are computed after excluding any NaN values. Y=MINMAX(X,N) with N between 0 and 1, returns generalized median filter value of X(:), e.g., it sorts elements of X then interpolates at (100*N) % rank position. N can be scalar, vector or matrix, so Y will have the same size as N. Examples: MINMAX(X,0) returns the minimum value of X elements. MINMAX(X,1) returns the maximum value of X elements. MINMAX(X,0.5) returns the median value of X elements. This is the equivalent of MEDIAN(X(:)). MINMAX(X,0.9) returns maximum value of X after excluding the 10% highest value elements. MINMAX(X,[0 1]) is the same as MINMAX(X). MINMAX(X,[0.01 0.99]) is a convenient way to compute automatic scale of X when samples are noisy, since it filters the 1% elements with extreme values. It may be used for color scaling with CAXIS or for plot scaling with set(gca,'YLim',...). See also MIN, MAX and MEDIAN. Author: François Beauducel <beauducel@ipgp.fr> Created: 1996 Updated: 20121010 Download minmax.m at Matlab Central file exchange (Copyright 2012 F. Beauducel / BSD License). 
This is the core formula of discrete Fourier transform: it simply computes the amplitude and phase shift of fundamental or harmonics of a phase signal. Example: t = linspace(0,2*pi); x = 2*cos(t + pi/2)  cos(3*t) + rand(size(t)); harmfit(t,x,1:4)returns estimations of amplitudes/phases for the first four harmonics. HARMFIT Sinusoidal harmonic curve fitting. H = HARMFIT(X,Y) returns first harmonic amplitude and phase of the data vector Y relative to phase vector X (in radian) in a 3element vector H = [N,AMP,PHA], so that AMP*COS(N*X + PHA) fits Y. HARMFIT(X,Y,N) computes the N'th harmonic. N = 1 stands for fundamental, N = 2 is second harmonic, etc... N can be a scalar or a vector of positive integers. For example use N = 1:4 to compute the first four harmonics. [H,YY] = HARMFIT(...) returns also an evaluation of harmonic curve fit in vector YY (as the sum of harmonics in N). This is simply the core calculation of discrete Fourier transform. Example: t = linspace(0,2*pi); x = 2*cos(t + pi/2)  cos(3*t) + rand(size(t)); harmfit(t,x,1:4) returns estimations of amplitudes/phases for the first four harmonics. Note that negative amplitudes are fitted with positive value with a PI phase difference. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20140522 Updated: 20140524 Download harmfit.m at Matlab Central file exchange (Copyright 2014 F. Beauducel / BSD License). 
These two scripts convert Roman numerals to and from any integer (scalar, vector or matrix), including large numbers greater than 4999 with the parenthesis notation (multiplies by 1000). The function NUM2ROMAN uses strict rules of modern notation (substractive principle for 4 and 9 bases) except for the common 'MMMM' form replacing '(IV)'. ROMAN2NUM is more flexible and is able to convert some other Roman notation possibilities, for instance the 3 different expressions roman2num({'IC','XCIX','XCVIIII'}) return [99,99,99], or roman2num({'MDXV','MCCCCCXV'}) return [1515,1515]. See help for syntax, and script comments for details. NUM2ROMAN Roman numerals. NUM2ROMAN(N) returns modern Roman numeral form of integer N, which can be scalar (returns a string), vector or matrix (returns a cell array of strings, same size as N). The function uses strict rules whith substractive notation and commonly found 'MMMM' form for 4000. It includes also parenthesis notation for large numbers (multiplication by 1000). Examples: num2roman(1968) num2roman(10.^(0:7)) reshape(num2roman(1:100),10,10) See also ROMAN2NUM. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 2005 Modified: 20091221 ROMAN2NUM Roman numerals conversion. ROMAN2NUM(X) converts to scalar the string X (or cell array of strings) of modern Roman numerals. The function considers strict rules with parenthesis notation for numbers larger than 4999. Invalid characters will produce NaN. Example: num2roman('MMMMDCCCLXXXVIII') num2roman('(CCCXIV)CLIX') See also NUM2ROMAN. Author: François Beauducel <beauducel@ipgp.fr> Institut de Physique du Globe de Paris Created: 20091219 Modified: 20091221 Download num2roman.m and roman2num.m at Matlab Central file exchange (Copyright 20052009 F. Beauducel / BSD License). 
This little function computes the date of Easter Sunday, using Oudin's algorithm for Julian calendar (starting 325 AD) and Gregorian calendar (after 1583 AD). Easter day is usefull to calculate other christian feasts, like Ash Wednesday (Easter  46), Good Friday (Easter  2), Ascension Thursday (Easter + 39), Pentecost (Easter + 49). Some translations of Easter: Pâques (French), Pascua (Spanish), Ostern (German), Pasqua (Italian), Paskah (Indonesian), Wielkanoc (Polish) ... See help for syntax, and script comments for details. EASTER Easter Day EASTER displays the date of Easter Sunday for present year. EASTER(YEAR) displays the date of Easter for specific YEAR, which can be scalar or vector. E = EASTER(...) returns Easter day(s) in DATENUM format. This function computes Easter Day using the Oudin's algorithm [1940], which is valid for Catholic Easter day from 325 AD (beginning of the Julian calendar). Easter day is usefull to calculate other usual christian feasts: datestr(easter47) is Mardi Gras datestr(easter46) is Ash Wednesday datestr(easter24) is MiCarême datestr(easter2) is Good Friday datestr(easter+1) is Easter Monday datestr(easter+39) is Ascension Thursday datestr(easter+49) is Pentecost Reference: Oudin, 1940. Explanatory Supplement to the Astronomical Almanac, P. Kenneth Seidelmann, editor. Tondering, C, 2008. http://www.tondering.dk/claus/calendar.html Author: François Beauducel, <beauducel@ipgp.fr> Created: 20021226 Updated: 20090301 Download easter.m at Matlab Central file exchange (Copyright 20022009 F. Beauducel / BSD License). 
The current Javanese calendar was inaugurated by Sultan Agung of Mataram in the Gregorian year 1633. It is based on a combination of the Hindu calendar "Saka" and the Islamic calendar based on the lunar month, and contains different cycles: Pasaran (5day), Dina Pitu (7day), Wetonan (35day), Mangsa (solar month), Wulan (Moon month), Pawukon (210day), Tahun (Moon year), Windu (8year), Kurup (120year). Coincidences of these multiple cycles have special mystical meanings for any Javanese people, for instance the birthday "Weton" or the Noble Days "Dino Mulyo". This is the primary timekeeping system for all matters having cultural, historical, and metaphysical significance in the Java island, Indonesia. This little script computes dates in the Javanese calendar, indicating Dina Pitu, Pasaran, Dino, Wulan, Tahun, Windu, and Kurup for today or a specific list of days. It indicates also the Noble Days if necessary. If you specify your date of birth, it can calculate your "weton" and a list of your javanese birthdays. For example, weton(1968,12,3) returns: A second use of this function is to display a month calendar that combines the 5day "Pasaran" cycle and the 7day Gregorian/Islamic week, called "Wetonan". See help for syntax, and script comments for details. WETON Javanese calendar / Wetonan. WETON without input argument returns the javanese date for today, in the form: WETON DINA WULAN T TAUN WINDU KURUP, DAY MONTH YEAR (DINO MULYO) where: WETON = DINAPITU PASARAN = combination of day names in the Gregorian/Islamic 7day week and Javanese 5day week D = day number in the Javanese month (1 to 29 or 30) WULAN = Javanese month name T = Javanese year number (starts on 1555) TAUN = Javanese year name (8 different, 12Wulan cycle) WINDU = Javanese "decade" name (4 different, 8Taun cycle) KURUP = Javanese "century" name (7 different, 120Taun cycle) DINO MULYO = "noble day" name (if necessary) WETON(YEAR,MONTH,DAY) returns the javanese date corresponding to YEARMONTHDAY in the Gregorian calendar. WETON(YEAR,MONTH,DAY,N) returns the list of your N first javanese birthdays (from the 35day "Weton" cycle). Example: if you are born on Dec 3, 1968 then weton(1968,12,3,10) returns your 10 first Wetons. Thanks to the Matlab flexibility, weton(1968,12,3+35*(0:10)) will do the same job... WETON(T) uses date T which can be Matlab date scalar, vector, matrix or any valid string for DATENUM function. Examples: weton(719135) weton('3Dec1968') weton([1968,12,3]) all return the string 'Selasa Kliwon 12 Pasa 1900 Ehé Adi Arbangiyah, 3 Desember 1968 (HANGGARA ASIH)' WETON(YEAR,MONTH) returns a javanese calendar for YEAR and MONTH in a table combining the 5day "Pasaran" cycle and 7day Gregorian week. Example: weton(1994,4) returns the following:  WETONAN BULAN APRIL 1994  Awal: Jemuwah Kliwon 19 Sawal 1926 Jé Sancaya Salasiyah, 1 April 1994 Akhir: Setu Wagé Dulkangidah 1926 Jé Sancaya Salasiyah, 30 April 1994  Senèn Selasa Rebo Kemis Jemuwah Setu Minggu Pon 04 19  14 29 09 24 Wagé 25 05 20  15 30 10 Kliwon 11 26 06 21 01 16  Legi  12 27 07 22 02 17 Pahing 18  13 28 08 23 03 where "Awal:" is the first day of the month, "Akhir:" the last one. X = WETON(...) returns a structure X with corresponding fields instead of displaying strings. To display full dates, use cat(1,char(X.weton)) Author: Mas François Beauducel Created: 19990127 (Rebo Pahing) Updated: 20110326 (Setu Pon) Download weton.m at Matlab Central file exchange (Copyright 19992011 F. Beauducel / BSD License). See also my little program weton.zip (Perl) that does almost the same job (and more). 